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Abstract 

The  purpose  of  this  thesis  is  to  examine  what  data  requirements  are  necessary 
to  avoid  continual  series  of  uncorrelated  tracks  when  gathering  observations.  The 
constants  of  the  motion  for  simple  two-body  motion  for  a  satellite  orbiting  the  Earth, 
known  as  the  classical  orbital  elements  or  COEs,  do  not  remain  constant  due  to 
zonal  and  sectoral  harmonic  variations  in  the  Earth’s  gravitational  field.  There  are 
other  elements  of  the  motion  that  should  be  considered  and  this  paper  discusses  the 
constancy  of  three  elements:  the  Hamiltonian  (TL)  of  the  Earth-Centered  Rotating 
System,  Z-component  of  inertial  angular  momentum  ( H A,  and  the  time  rate  of  change 
of  the  right  ascension  of  the  ascending  node  (f2). 

With  an  understanding  of  the  constancy  of  these  elements,  simulated  data  was 
used  to  determine  the  effects  sensor  performance  and  observation  quantity  have  on  the 
ability  to  effectively  estimate  these  constants.  This  information  was  used  to  determine 
an  appropriate  level  of  fidelity  for  a  model  to  be  utilized  as  a  supplement  in  fitting 
observation  data  with  current  data  available  in  the  Satellite  Catalog. 
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UNCORRELATED 
TRACK  AVOIDANCE 

I.  Introduction 

1.1  Background, 

Since  the  first  man-made  satellite  Sputnik  was  launched  by  the  former  Soviet 
Union  in  1957,  the  number  of  Earth  orbiting  satellites  continues  to  grow  at  a  rapid 
pace.  Currently,  the  Earth  is  surrounded  by  thousands  of  satellites  and  other  objects, 
such  as  expended  rocket  bodies  and  debris  left  from  broken-up  satellites  orbiting  at 
numerous  altitudes  above  the  Earth.  These  objects  create  a  region  of  increased  risk 
around  the  Earth  for  objects  returning  into  the  Earth’s  atmosphere  or  objects  being 
launched  into  space. 

These  problems  are  not  new  to  anyone  in  the  space  business  but  they  do  re¬ 
quire  full  time  operations  to  monitor  and  track  all  of  these  objects.  Currently,  the 
United  States  uses  the  “Satellite  Catalog  (SATCAT)”  [6]  to  track  over  8,500  objects 
approximately  as  small  as  a  softball  [1],  In  order  to  track  these  different  objects, 
current  operations  entail  obtaining  observational  data  and  comparing  the  new  data 
with  databases  filled  with  archived  information  on  current  objects.  Then,  based  on 
best-fit  algorithms  like  least  squares,  the  tracked  satellite  is  fit  to  an  existing  track  in 
the  database  and  the  information  for  that  object  is  updated;  if  the  new  orbit  cannot 
be  fit  with  an  existing  object  in  the  database,  it  is  labeled  as  an  uncorrelated  track 
and  inserted  into  the  database  as  a  new  object.  In  the  business  where  accuracy  is  of 
utmost  importance  and  information  is  essential,  it  is  beneficial  to  keep  uncorrelated 
tracks  at  a  minimum.  Tracking  all  of  these  objects  is  quite  a  task  especially  with 
all  of  the  disturbing  forces  that  are  acting  on  these  objects,  varying  with  each  orbit. 
Despite  the  difficulty  of  the  task  and  the  standing  army  required  to  manage  the  prob¬ 
lem,  there  are  ideas  to  reduce  the  size  of  the  objects  being  tracked  from  softball  size 


1 


to  approximately  golf  ball  size  which  will  substantially  increase  the  number  of  objects 
being  tracked. 

1.2  Problem  Statement 

As  the  number  of  objects  being  tracked  increases,  one  can  naturally  expect 
that  the  number  of  uncorrelated  tracks  will  increase  at  the  same  rate.  However, 
the  problem  is  not  based  on  numbers  alone  and  because  the  objects  being  tracked 
are  smaller  in  size,  the  objects  may  only  be  in-view  when  they  pass  overhead  or 
nearly-overhead.  Based  on  the  orbit  of  the  object  and  the  location  of  the  tracking 
station,  these  objects  might  only  be  tracked  once  every  two  weeks.  That  means  the 
information  in  the  catalog  may  not  be  current  when  the  object  is  viewed  again,  giving 
rise  for  more  error  in  the  attempt  to  fit  the  new  track  with  one  in  the  database.  The 
higher  error  increases  the  likelihood  of  adding  identical  objects  multiple  times  into 
the  catalog.  Therefore,  in  order  to  limit  the  redundancy  in  the  Satellite  Catalog  by 
minimizing  the  number  of  uncorrelated  tracks,  there  is  a  need  to  enhance  the  process 
by  including  information  about  the  orbit  that  remains  more  accurate  for  a  longer 
period  of  time. 

1.3  Research  Focus 

The  focus  of  this  research  is  to  investigate  the  constancy  of  the  Hamiltonian 
Function  (Chapter  III),  the  constancy  of  the  Z-component  of  the  inertial  angular 
momentum  (Chapter  IV)  and  the  constancy  of  the  time  rate  of  change  of  the  right 
ascension  of  the  ascending  node  (Chapter  V).  The  constancy  of  these  elements  is 
established  through  analysis  of  different  levels  of  fidelity  for  the  gravitational  geopo¬ 
tential  expansion.  The  level  of  accuracy  attainable  for  each  of  these  elements  is  then 
utilized  in  establishing  a  basis  for  a  supplemental  model  (Chapter  VI)  that  can  be 
used  in  data  mining  to  aide  in  the  minimization  of  uncorrelated  tracks  by  fitting  two 
or  more  uncorrelated  tracks. 
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II.  Coordinate  Systems 


2.1  Earth- Centered  Coordinate  Frames 

There  are  two  coordinate  systems  that  are  referenced  in  this  paper  and  a  short 
description  of  each  coordinate  system  is  addressed  in  this  chapter.  The  first  coordinate 
frame  is  the  Earth-Centered  Inertial  (ECI)  Frame,  also  known  as  the  Geocentric 
Equatorial  Coordinate  System  or  UK  Frame.  The  origin,  as  suggested  in  the  name, 
lies  at  the  center  of  the  Earth.  The  fundamental  plane  for  this  system  is  the  Earth’s 
equatorial  plane  and  the  first  axis  (/)  points  toward  the  vernal  equinox  or  the  first 
point  of  Aries.  The  second  axis  ( J)  is  picked  ninety  degrees  off  of  the  first  and  chosen 
so  that  applying  the  right-hand  rule,  the  third  component  of  the  coordinate  system 
(K)  is  aligned  with  the  Earth  angular  velocity  vector,  through  the  North  Pole  [8,  157]. 
Figure  2.1  shows  the  ECI  Coordinate  Frame  and  comes  from  Vallado  [8,  157]. 


Figure  2.1:  Earth-Centered  Inertial  Frame 


The  second  coordinate  system  that  is  referenced  in  this  paper  is  the  Earth- 
Centered  Rotating  (ECR)  Frame,  also  known  as  the  Body-Fixed  Coordinate  System 
or  International  Terrestrial  Reference  Frame  (ITRF).  This  coordinate  system  shares 
the  origin  and  fundamental  plane  of  the  ECI,  but  the  main  difference  is  that  the 
coordinate  system  rotates  with  the  Earth.  The  Erst  axis  (A")  points  from  the  origin 
toward  the  Prime  Meridian,  or  Greenwich  Meridian,  which  is  zero  degrees  longitude. 
Similar  to  the  ECI,  the  second  axis  (E)  is  selected  ninety  degrees  East  of  first  axis 
such  that  the  third  axis  (Z)  is  aligned  with  the  Earth  spin  axis  using  the  right-hand 
rule. 


3 


As  stated  above,  one  difference  between  these  two  coordinate  systems  is  that 
ECR  rotates  with  the  Earth  and  a  link  between  these  two  coordinate  frames  is  the 
angle  known  as  Greenwich  Meridian  Sidereal  Time  ( Oqmst )■  Greenwich  Meridian 
Sidereal  Time  is  the  angle  measured  from  the  first  axis  of  the  ECI  (/)  to  the  Earth’s 
meridian  [5,  137].  Figure  2.2  is  taken  from  Vallado  [8,  189]  but  some  of  the  information 
was  removed  in  order  to  focus  on  the  angles  of  interest.  Additionally,  due  to  slight 


Figure  2.2:  Greenwich  Meridian  Sidereal  Time 

movements  over  time  by  the  equatorial  plane  and  the  equinox,  nutation  and  precession 
models  are  also  required  to  link  the  two  coordinate  systems  [8,  158]. 
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III.  Integrals  of  Satellite  Motion 


3. 1  Introduction 

The  first  element  the  author  analyzed  for  constancy  is  the  Hamiltonian  Func¬ 
tion.  Solving  for  the  Hamiltonian  Function  is  accomplished  by  utilizing  analytical 
mechanics.  The  advantages  to  using  analytical  mechanics  is  that  instead  of  focus¬ 
ing  on  vector  quantities  like  force  and  momentum,  the  focus  shifts  to  scalars  like 
kinetic  energy  and  the  Hamiltonian.  Not  only  does  changing  from  vectors  to  scalars 
make  things  easier,  but  the  acceleration  vector  is  not  required  in  analytical  mechanics. 
Additionally,  analytical  mechanics  allows  the  use  of  generalized  coordinates  allowing 
selection  of  an  ideal  coordinate  system  to  calculate  the  above  mentioned  scalars.  [7,  45] 


3.2  Deriving  the  Hamiltonian  Function 

In  formulating  the  Hamiltonian  Function,  x,  y,  and  z  in  the  Earth-Centered 
Rotating  Frame  are  chosen  to  be  the  “generalized  coordinates”  ( q )  [7,  45]: 


ql 

X 

Q  = 

<?2 

= 

y 

V  93  ) 

V  *  / 

(3.1) 


The  next  step  is  to  take  the  derivative  of  the  generalized  coordinates  to  find  the 
“generalized  velocities”  [7,  94]  and  because  the  q'ks  are  in  a  rotating  frame,  we  must 
apply  the  following  equation,  taken  from  Wiesel  [9,  12]  to  find  the  inertial  derivative: 
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Note  that  in  the  above  equation  usl  is  ca  of  the  Earth  and  it  is  labeled  as  cu®.  Taking 
the  derivative  of  the  coordinates  produces  the  following  generalized  velocities  (q): 


f9Il 
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= 

y 

+  x 

92 

U3 ) 

u  > 
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Executing  the  cross  product  in  the  above  equation  yields  the  following: 
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Now  obtaining  the  generalized  velocities  makes  it  possible  to  calculate  the  kinetic 
energy  (T)  by  using  the  following  equation  [7,  14]: 


f91l 
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(3.5) 


However,  because  the  interest  is  in  the  specific  kinetic  energy,  there  is  no  need  for 
mass  in  the  above  equation.  Taking  out  the  mass  term  and  inputting  the  results  from 
Equation  3.4  the  specific  kinetic  energy  is  calculated  to  be: 


T 


1 

2 


(x  -  uaq2f  +  (jl  +  u>®9i  )2  +  (zf 


(3.6) 


Along  with  the  specific  kinetic  energy  of  the  system,  the  potential  energy  (V)  of  the 
system  is  also  of  interest.  The  potential  energy  can  be  written  in  the  geopotential 
form  with  r  =  y/gf-|-g|~Tg|  as  follows  [10,  108]: 


v  =  ~^yy 

nr  /  J  /  J 


n= 0  m=0 


r 


p:n  (—  ix 


Cnrn  cos  m  tan 


Qi 


.  -1/92 

nm  svri  m  tan  — 

(1\ 


(3.7) 
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The  Cnm  and  Snm  terms  were  defined  in  a  study  conducted  by  a  team  comprised  of  in¬ 
dividuals  from  the  National  Geospatial-Intelligence  Agency  and  National  Aeronautics 
and  Space  Administration.  The  study  assigned  values  for  the  geopotential  expansion 
out  through  n  —  m  —  50  [2],  The  final  n  and  m  values  desired  for  the  geopotential 
expansion  will  be  referred  to  as  J™  for  the  remainder  of  this  paper.  This  provides  the 
n  and  m  values  required  for  the  summations  to  determine  the  values  for  Cnm,  Snm 
and  P™  in  Equation  3.7. 

After  obtaining  values  for  the  specific  kinetic  energy  and  potential  energy  of  the 
system,  the  next  step  is  to  find  the  Lagrangian  (L)  [7,  68].  The  Lagrangian  is  defined 
as: 

L  =  T  —  V  (3.8) 

Inserting  Equation  3.6  and  3.7  into  Equation  3.8  the  Lagrangian  is  found  to  be: 

L  =  \  [0*  j  u®q2f  +  (y  +  w©gi)2  +  (i)2]  -  V (qi,  q2,  q3)  (3.9) 


The  next  step  is  to  calculate  the  “generalized  momenta”  (p)  [7,  80].  The  momenta 
are  defined  as: 


dL  dT 
Pk  —  -7TT-  ~  WW 
OQk  Oqk 

Applying  equation  3.10  it  is  possible  to  find  the  generalized  momenta: 


(3.10) 


/  \ 

(  '  \ 

Pi  ' 

x  -  cneg2 

P2 

= 

V  +  ueqi 

\P3  ) 

V  i  / 

(3.11) 


Which  shows  that  the  generalized  momenta  are  the  same  as  the  inertial  derivative  of 
the  coordinates  or  the  generalized  velocities.  Therefore,  the  momenta  for  this  non- 
inertial  system  are  the  components  of  the  inertial  velocity.  Rearranging  the  above 
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equation  to  solve  for  the  qk's  the  following  is  found: 


/  .  \ 

/  ■  \ 

(  , 

\ 

<1\  ’ 

X 

Pi  +  ueq2 

4-2 

= 

y 

— 

P2  —  w ©<7i 

{43 ) 

V  i  ) 

\  P3 

) 

(3.12) 


The  last  step  is  to  formulate  the  Hamiltonian  [7,  94],  The  Hamiltonian  Function  is 
defined  as: 

n 

H  =  ^pkqk~L  (3.13) 

k=  1 

but  the  generalized  velocities  must  not  be  in  the  Hamiltonian  Function  [7,  94];  there¬ 
fore,  substituting  in  for  the  generalized  velocities,  as  defined  in  Equation  3.12,  results 
in  a  Hamiltonian  Function  as  defined  below: 


n 


1 

2 


(pi  +  vl  +  Pa)  +  (Pi«2  -  P29i)  +  V (ji ,  <)2, 93) 


(3.14) 


Equation  3.14  shows  that  the  Hamiltonian  Function  is  not  directly  dependent  on  time 
making  it  a  constant  of  the  motion. 

It  is  worth  noting  that  the  above  expression  can  be  simplified  by  incorporating 
polar  coordinates.  In  polar  coordinates,  r  is  the  radius  as  defined  above,  4>  is  defined  as 
the  geocentric  longitude  and  9  is  defined  as  the  colatitude  [10,  105].  These  coordinates 
are  related  to  the  generalized  coordinates  in  the  following  manner: 

<t>  =  tan,-1  ^  (3.15) 


9  =  cos  1 


93 


(3.16) 


Substituting  Equations  3.15  and  3.16  into  Equation  3.7  results  in  the  Hamiltonian 
Function  expressed  in  polar  coordinates  as  defined  in  Wiesel  [10,  108].  Now  working 
with  the  Hamiltonian  Function  in  polar  coordinates  it  is  possible  to  calculate  the 


change  in  the  generalized  momenta  using  the  following  equation  [7,  94]: 


dn 

Pk  =  -x— 

oqk 


(3.17) 


Taking  the  partial  as  described  in  3.17  results  in  a  change  in  the  generalized  momenta 
for  the  geocentric  longitude  (0)  as  shown  below  [10,  110]: 

oo  n  /  \  —  n 

P4>  =  ~  V!  V!  f  TT  )  pn  ( cos6)xm  [—Cnm sinm  (0  -  u®t)  +  Snmcosm  (0  -  cuet)] 

r  n=0  m=0 

(3.18) 

This  result  shows  that  for  cases  where  m  —  0  the  change  in  p $  is  equal  to  zero 
making  the  inertial  angular  momentum  a  constant  of  the  motion.  This  is  the  topic  of 
discussion  for  Chapter  IV. 


3.3  Constancy  of  the  Hamiltonian 

In  order  to  show  the  constancy  of  the  Hamiltonian  Function,  the  author  used 
Matlab®  to  create  a  computer  program  (see  Appendix  B)  which  takes  position  and 
velocity  vectors  in  the  Earth-Centered  Rotating  (ECR)  Frame  and  computes  the 
Hamiltonian  as  described  in  Equation  3.14  using  simulated  data.  The  following  is  a 
complete  list  of  the  tests  run  to  determine  the  constancy  of  the  Hamiltonian  Function: 


•  Case  1: 

•  Case  2: 

•  Case  3: 

•  Case  4: 

•  Case  5: 

•  Case  6: 

•  Case  7: 

•  Case  8: 


Includes  geopotential  expansion  through  term 
Includes  geopotential  expansion  through  term 
Includes  geopotential  expansion  through  J|  term 
Includes  geopotential  expansion  through  term 
Includes  geopotential  expansion  through  Jj  term 
Includes  geopotential  expansion  through  term 
Includes  geopotential  expansion  through  Jf  term 
Includes  geopotential  expansion  through  J9n  term 
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•  Case  9:  Includes  geopotential  expansion  through  J^1  term 

The  above  cases  were  run  using  four  different  scenarios.  The  first  scenario  is  a 
near-circular  orbit,  inclined  at  30  degrees,  evaluated  for  2000  TU  or  approximately 
18  days.  The  second  scenario  is  a  circular  orbit,  inclined  at  45  degrees,  evaluated  for 
2000  TU.  The  third  scenario  is  the  same  orbit  as  Scenario  1  with  an  evaluation  period 
of  20000  TU  or  approximately  six  months.  The  final  scenario  is  the  same  orbit  as 
Scenario  2  with  an  evaluation  period  of  20000  TU.  The  initial  conditions  for  position 
( Rijk )  and  velocity  ( Vjjk )  f°r  Scenarios  1  and  3  are: 

1.05 
0.0 
0.0 

0.0 

Vjjk  =  0.845154254 

v  0.453464763 

The  initial  conditions  for  Scenarios  2  and  4  are: 

(  0.0 

Rijk  =  -1.15678606465 

0.0 

0.65744356375 

0.0 

0.65744356375 

For  each  of  the  nine  cases  and  four  scenarios  the  process  of  gathering  the  data 
is  the  same.  The  desired  initial  conditions  and  geopotential  value  are  input  into  Dr. 
Wiesel’s  orbit  propagator  by  altering  the  input  hie  to  specify  the  initial  conditions 
as  well  as  the  n  and  m  values  for  the  geopotential  expansion.  The  orbit  propagator 
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outputs  position  and  velocity  data  in  both  the  ECI  and  ECR  Frames  which  are  then 
read  into  the  Matlab®  program.  Then  data  showing  all  the  values  for  the  Hamiltonian 
for  each  term  at  each  time  step  is  output. 

Applying  the  initial  conditions  for  Scenario  1,  the  first  test  case  was  evaluated 
in  which  n  =  m  =  0  to  find  the  two-body  constancy  of  the  Hamiltonian  Function. 
Taking  the  output  data  from  the  computer  program  and  importing  it  into  Excel® 
a  graph  of  the  values  at  each  time  step  was  created.  Because  this  case  is  the  two- 
body  case,  there  are  no  geopotential  terms  to  alter  the  constancy  of  the  system  and, 
therefore,  a  constant  value  for  H  for  the  entire  data  collection  period  is  expected. 


Hamiltonian  Function  -  J0,0 


Figure  3.1:  Scenario  1,  Case  1:  H  versus  time 

Figure  3.1  verifies  that  indeed  the  Hamiltonian  Function  is  constant  up  to  ap¬ 
proximately  order  10~13.  Because  the  value  of  the  Hamiltonian  Function  is  calculated 
using  canonical  units,  the  original  values  are  of  order  one,  making  the  7i  constant  to 
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twelve  significant  figures.  This  result  was  expected  for  the  two-body  case  and  will  be 
used  as  a  comparison  for  the  remainder  of  the  cases. 

The  next  case  evaluated  was  Case  2  with  geopotential  terms  out  to  J°.  Using 
the  same  process  stated  above  data  was  gathered  and  the  results  are  plotted  in  Figure 
3.2.  This  figure  shows  that  the  constancy  of  7i  remains  at  approximately  10~13  even 
with  the  oblateness  of  the  Earth  factored  into  the  equation.  The  importance  of  these 
results  is  that  although  the  geopotential  term  is  added  the  constancy  of  H  remains 
on  the  same  order  as  the  two-body  case.  When  analyzing  the  COEs,  there  is  a 
noticeable  change  in  the  accuracy  of  these  constants  with  the  oblateness  term  added 
and  equations  for  the  secular  change  can  be  found  in  Wiesel  [10,  141-146].  Unlike 
the  classical  orbital  elements  though,  TL  retains  the  same  order  of  constancy  with  the 
inclusion  of  the  gravity  geopotential  expansion  through  the  J3  term. 


Hamiltonian  Function  -  J2,0 


Figure  3.2:  Scenario  1,  Case  2:  7i  versus  time 
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Continuing  to  add  terms  from  the  geopotential  does  not  affect  the  accuracy 
of  the  Hamiltonian  Function.  Executing  all  cases  in  the  Erst  scenario  results  in  no 
change  to  the  level  of  constancy  of  Hi.  The  constancy  of  7i  remains  at  approximately 
order  ICC13.  Scenario  1,  Case  9  validates  that  adding  geopotential  terms  through  the 
J2I  term  the  Hamiltonian  Function  remains  at  a  constancy  of  approximately  the  same 
order  of  magnitude  as  the  two-body  case.  Figure  3.3  is  a  graph  of  Scenario  1,  Case 
9  graphically  showing  that  for  2000  TU,  7i  remains  constant  to  approximately  order 
ICC13.  As  was  stated  above  all  of  the  cases  in  this  scenario  have  approximately  the 
same  constancy  as  Case  1  and  the  graphs  for  Scenario  1  not  explicitly  discussed  in 
this  section  are  attached  in  Appendix  A  along  with  graphs  for  Scenario  2. 


Hamiltonian  Function  -  J21.21 


Figure  3.3:  Scenario  1,  Case  9:  H  versus  time 

In  addition  to  graphing  the  data,  analysis  was  done  in  order  to  End  the  maximum 
change  from  the  initial  Hamiltonian  Function  value.  In  order  to  do  this,  the  data  was 
inserted  into  an  Excel®  file  and  the  absolute  value  of  the  difierence  of  lH.it)  —  7d(io) 
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was  calculated  in  one  column.  Then  using  a  built-in  Excel®  function  the  maximum 
value  in  that  column  was  identified.  For  each  of  the  four  scenarios  and  nine  cases  in 
each  scenario,  the  data  was  evaluated  in  the  same  manner  and  the  results  from  all  of 
these  cases  are  summarized  in  Table  3.1.  Note  the  asterisk  next  to  Scenario  3,  Case  1; 
it  is  there  to  signify  that  data  was  not  collected  for  the  entire  20000  TU.  For  unknown 
reasons,  the  orbit  propagator  did  not  produce  data  for  the  entire  time  interval  and 
therefore  the  data  is  based  on  a  shorter  time  period.  This  is  insignificant  since  it  is 
only  the  two-body  case  and  we  expect  it  to  remain  constant  for  the  entire  time. 


Scenario  1 

Scenario  2 

Scenario  3 

Scenario  4 

Case  1 

2.46E-13 

8.80E-14 

3.29E-13* 

6.29E-13 

Case  2 

2.18E-13 

4.68E-13 

1.51E-12 

9.73E-13 

Case  3 

4.56E-13 

1.13E-13 

8.19E-13 

5.79E-13 

Case  4 

4.17E-13 

2.44E-13 

3.96E-13 

9.52E-13 

Case  5 

2.92E-13 

1.15E-13 

7.99E-13 

1.10E-12 

Case  6 

4.16E-13 

8.20E-14 

1.44E-12 

7.37E-13 

Case  7 

1.79E-13 

2.15E-13 

6.67E-12 

5.45E-13 

Case  8 

2.86E-13 

1.23E-13 

1.35E-12 

3.11E-13 

Case  9 

1.35E-13 

1.68E-13 

4.18E-13 

4.39E-13 

Table  3.1:  Hamiltonian  Function  Constancy  Results 


The  final  analysis  done  with  the  Hamiltonian  Function  was  to  come  up  with  a 
comparison  for  how  well  a  lower  geopotential  expansion  estimate  does  with  a  specific 
data  set.  The  motivation  behind  this  analysis  is  to  determine  a  level  of  constancy 
for  the  Hamiltonian  Function  based  on  position  and  velocity  accuracy  attainable 
for  a  specific  sensor.  For  example,  it  may  not  be  feasible  to  use  the  geopotential 
expansion  through  the  J term  if  the  position  accuracy  for  a  sensor  is  on  the  order 
of  a  kilometer.  This  analysis  will  provide  a  basis  for  the  discussion  in  Chapter  VI  to 
determine  an  appropriate  geopotential  expansion  level  for  models  based  on  position 
and  velocity  accuracy  attainable  by  a  sensor.  In  order  to  accomplish  this,  position 
and  velocity  data  was  created  by  running  Dr.  Wiesel’s  orbit  propagator  for  each  of 
the  four  scenarios  detailed  above.  The  geopotential  input  into  the  propagator  calls  for 
the  n  =  21  and  m  =  21  values  indicating  the  geopotential  expansion  through  the  Jf* 
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term.  As  discussed  previously,  the  output  from  the  orbit  propagator  is  used  as  input 
into  the  Hamiltonian  Function  Matlab®  code.  Then  Hamiltonian  data  is  gathered 
to  show  the  accuracy  of  each  of  the  different  cases  as  defined  above.  So  for  Case  1, 
the  Hamiltonian  data  created  for  the  geopotential  expansion  through  the  J|4  term 
was  evaluated  using  only  the  two-body  terms  in  the  Matlab®  program.  Then  taking 
the  output  from  the  computer  program  and  importing  it  into  Excel®  the  data  was 
analyzed  as  described  above.  A  column  of  data  was  created  containing  the  absolute 
value  of  the  difference  between  H(t)  —  7i(t0)  and  then  using  the  built-in  Excel® 
function  the  maximum  difference  from  the  initial  Hamiltonian  value  was  calculated 
and  inserted  into  Table  3.2.  The  accuracy  of  the  Hamiltonian  Function  with  Case  1 
assumptions  was  very  low,  on  the  order  of  10”4.  The  case  increases  the  constancy 
of  7i  by  about  an  order  of  magnitude,  but  is  still  not  very  constant.  It  turns  out  that 
Case  2  is  about  as  accurate  as  the  “zonal  harmonic”  cases  will  get  [10,  113].  Note  in 
the  table  that  using  ./^  yields  approximately  the  same  constancy.  However,  it  should 
be  noted  that  as  more  terms  are  added  to  the  geopotential  in  the  Matlab®  program, 
the  constancy  of  the  Hamiltonian  Function  increases.  These  results  are  summarized 
in  Table  3.2  and  the  table  ends  with  the  constancy  obtained  by  using  the  geopotential 
through  Jf4  exactly  what  was  used  to  create  the  position  and  velocity  vector  in  the 
orbit  propagator.  The  final  results  in  Table  3.2  match  the  constancy  in  the  final  row 
of  Table  3.1  because  the  data  input  into  the  propagator  and  Matlab®  program  is 
identical. 
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Scenario  1 

Scenario  2 

Scenario  3 

Scenario  4 

Case  1 

3.23E-04 

5.23E-04 

3.23E-04* 

5.25E-04 

Case  2 

2.12E-05 

1.23E-05 

2.13E-05 

1.24E-05 

Case  3 

1.15E-05 

7.28E-06 

1.24E-05 

7.30E-06 

Case  4 

2.11E-05 

1.23E-05 

2.11E-05 

1.23E-05 

Case  5 

8.17E-06 

2.56E-06 

8.20E-06 

2.58E-06 

Case  6 

2.12E-05 

1.23E-05 

2.12E-05 

1.23E-05 

Case  7 

3.46E-06 

6.03E-07 

3.93E-06 

6.03E-07 

Case  8 

2.13E-05 

1.23E-05 

2.13E-05 

1.23E-05 

Case  9 

1.35E-13 

1.68E-13 

4.18E-13 

4.39E-13 

Table  3.2:  Hamiltonian  Function  Comparison 
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IV.  Angular  Momentum 


4  ■  1  Introduction 

Another  component  evaluated  for  constancy  is  the  Z-component  of  inertial  an¬ 
gular  momentum  of  the  system.  Starting  in  the  Earth-Centered  Inertial  Frame  and 
neglecting  all  disturbing  forces  other  than  the  zonal  harmonics  in  the  geopotential, 
the  inertial  angular  momentum  stays  constant  with  a  symmetric  Earth,  this  result 
was  proved  in  Chapter  III.  Inertial  angular  momentum  ( Hijk )  is  defined  as: 


Hijk  =  Rijk  X  Vjjk  (4.1) 

As  noted  above,  as  long  as  the  Earth  is  symmetric  about  the  polar  axis,  the  inertial  an¬ 
gular  momentum  of  the  system  will  stay  constant.  The  constancy  of  the  Z-component 
of  inertial  angular  momentum  was  evaluated  using  a  Matlab®  program  (see  Appendix 
B).  This  program  reads  in  the  ECI  position  and  velocity  vectors  and  computes  the 
angular  momentum  at  each  time  step.  As  noted  in  the  previous  chapter,  the  position 
and  velocity  vectors  are  calculated  using  Dr.  Wiesel’s  custom-made  orbit  propagator. 

4-2  Changing  the  Zonal  Harmonics 

In  order  to  quantify  the  constancy  of  the  Z-component  of  inertial  angular  mo¬ 
mentum  analysis  of  the  same  four  scenarios  as  stated  in  Chapter  III  was  accomplished. 
The  four  cases  include  the  near-circular  orbit  with  a  thirty  degree  inclination  and  the 
circular  orbit  with  a  forty-five  degree  inclination  evaluated  for  2000  TU  and  20000 
TU  for  each  of  these  two  orbits.  Additionally,  each  of  these  four  scenarios  was  ana¬ 
lyzed  using  the  same  nine  cases  as  discussed  in  the  previous  chapter,  but  the  analysis 
was  done  in  a  different  order.  The  modification  to  the  order  is  used  to  point  out 
the  difference  between  a  symmetric  and  non-symmetric  Earth.  Initially,  analysis  for 
the  scenarios  in  which  only  the  “zonal  harmonics”  are  considered  was  accomplished, 
keeping  the  Earth’s  gravitational  geopotential  symmetric  [10,  113].  The  five  initial 
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cases  run  on  each  scenario  to  determine  the  constancy  of  the  Z-component  of  the 
inertial  angular  momentum  are  listed  below: 

•  Case  1:  Includes  geopotential  expansion  through  the  Jg  term 

•  Case  2:  Includes  geopotential  expansion  through  the  term 

•  Case  4:  Includes  geopotential  expansion  through  the  term 

•  Case  6:  Includes  geopotential  expansion  through  the  Jg  term 

•  Case  8:  Includes  geopotential  expansion  through  the  J9n  term 

Similarly,  the  initial  conditions  for  the  two  scenarios  remain  the  same  as  stated  in 
Chapter  III.  Initial  conditions  for  Scenarios  1  and  3  are: 

f  1.05 
Rijk  =  0.0 

v  0.0 

(  0.0 

Vijk  =  0.845154254 

v  0.453464763 

Initial  conditions  for  Scenarios  2  and  4  are: 

0.0 

-1.15678606465 

0.0 

0.65744356375 

0.0 

0.65744356375 

The  first  case  analyzed  is  Scenario  1,  Case  1,  where  n  —  m  —  0  in  the  geopo¬ 
tential  (Equation  3.7).  In  order  to  verify  the  constancy,  the  Z-component  of  the 


[DU]  (4,4) 
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inertial  angular  momentum,  (Hk),  is  plotted  against  time  represented  in  canonical 
time  units  [TU] .  This  was  accomplished  by  calculating  inertial  angular  momentum 
using  a  Matlab®  program  and  the  Z-component  values  were  output  along  with  each 
associated  time.  The  data  was  then  imported  into  Excel®  and  plotted,  showing  the 
constancy  of  Hk  versus  time. 


Specific  Angular  Momentum  -  J0,0 


Time  [TU] 


Figure  4.1:  Scenario  1,  Case  1:  Hk  versus  time 

Figure  4.1  shows  that  for  the  initial  conditions  stated  in  Equations  4.2  and 
4.3,  the  Z-component  of  the  inertial  angular  momentum  is  constant  to  approximately 
order  10-14  over  a  period  of  2000  TU,  or  just  over  two  weeks.  This  was  expected  since 
the  first  case  is  merely  the  two-body  case  without  the  effects  of  the  geopotential, 
the  Z-component  of  the  inertial  angular  momentum  is  constant  along  with  all  the 
other  classical  orbital  elements.  This  first  case  will  be  used  for  comparison  with  the 
remaining  cases  within  this  scenario. 
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Specific  Angular  Momentum  -  J2,0 


Figure  4.2:  Scenario  1,  Case  2:  Hk  versus  time 

The  next  case  evaluated  is  Case  2,  which  accounts  for  the  oblateness  of  the 
Earth,  or  the  Earth’s  “equatorial  bulge”  [10,  114].  Running  the  same  Matlab®  pro¬ 
gram  as  before  and  importing  the  data  into  Excel®,  a  similar  graph  is  created.  Figure 
4.2  plots  the  Z-component  of  the  inertial  angular  momentum  against  time  and  shows 
the  constancy  remaining  at  approximately  order  10~13  with  the  stated  initial  condi¬ 
tions.  This  continues  to  validate  the  assertion  that  the  Z-component  of  the  inertial 
angular  momentum  remains  constant  over  a  symmetric  Earth.  This  first  zonal  change 
in  the  geopotential  creates  a  symmetric  change  in  the  Earth’s  gravitational  held  there¬ 
fore,  leaving  the  Z-component  of  the  inertial  angular  momentum  constant. 

This  trend  continues  for  all  cases  where  the  geopotential  contains  no  “sectoral 
harmonics”  or  “tesseral  harmonics”,  or  where  m  =  0  [10,  114-115].  Figure  4.3  is  a 
graph  of  the  results  calculated  for  the  J3,  case  and  it  shows  that  the  constancy  of  the 
Z-component  of  the  inertial  angular  momentum  continues  to  remain  at  a  constant 
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Specific  Angular  Momentum  -  J21 ,0 


Figure  4.3:  Scenario  1,  Case  8:  Hk  versus  time 

value  to  within  approximately  order  10~13.  Case  8  is  the  final  case  that  was  ana¬ 
lyzed  to  verify  the  constancy  of  the  Z-component  of  the  inertial  angular  momentum 
for  Scenario  1.  The  constancy  of  this  case  is  unchanging  from  the  previous  cases 
including  the  two-body  case  showing  that  for  a  symmetric  Earth,  the  Z-component 
of  the  inertial  angular  momentum  will  remain  constant  to  approximately  order  10-13. 
Graphs  from  Scenario  2  and  cases  not  cited  in  this  section  can  be  found  in  Appendix 
A. 

4-3  Including  the  Tesseral  and  Sectoral  Harmonics 

Although  making  use  of  the  assumption  that  the  Earth’s  gravitational  field 
is  symmetric  allows  for  great  constancy  in  the  Z-component  of  the  inertial  angular 
momentum,  it  does  not  represent  the  Earth’s  actual  gravitational  field.  Including 
the  tesseral  and  sectoral  harmonics  in  the  calculations  enhances  our  ability  to  more 
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accurately  model  the  actual  effects  that  the  Earth’s  geopotential  has  on  orbiting 
objects.  Now  utilizing  the  same  four  scenarios  as  above,  the  constancy  of  the  Z- 
component  of  the  inertial  angular  momentum  was  evaluated  using  the  following  four 
cases: 

•  Case  3:  Includes  geopotential  expansion  through  the  Jf  term 

•  Case  5:  Includes  geopotential  expansion  through  the  J|  term 

•  Case  7:  Includes  geopotential  expansion  through  the  Jf  term 

•  Case  9:  Includes  geopotential  expansion  through  the  Jf]1  term 

Each  of  these  cases  includes  the  geopotential  expansion  through  a  term  where  n  =  m. 
These  terms  are  referred  to  as  “sectoral  harmonics”  [10,  114]. 

The  first  case  to  be  analyzed  is  Case  3.  In  Scenario  1,  the  near-circular,  30 
degree  inclined  orbit,  was  evaluated  by  adding  tesseral  and  sectoral  harmonics  up  to 
Jf .  A  graph  was  created  to  show  the  value  of  the  Z-component  of  the  inertial  angular 
momentum  versus  time.  The  data  collected  from  Matlab®  was  imported  into  Excel® 
and  plotted  for  2000  TU.  Figure  4.4  shows  that  the  constancy  of  the  Z-component 
of  the  inertial  angular  momentum  decreases  from  order  ICC14  as  seen  in  Case  1  and 
order  ICC13  as  seen  in  Case  8  to  a  constancy  on  the  order  of  approximately  ICC5. 
This  is  a  huge  reduction  in  constancy  and  is  due  solely  to  the  non-symmetric  nature 
of  the  tesseral  and  sectoral  harmonics.  It  can  be  noted  in  Figure  4.4  that  each  time 
the  object  orbits  the  Earth  there  are  periodic  effects  causing  the  Z-component  of  the 
inertial  angular  momentum  to  vary  equally  over  the  same  location  of  the  Earth.  It 
appears  that  the  period  of  these  cyclic  variances  is  related  to  the  period  of  the  orbit. 

The  next  case  evaluated  in  this  scenario  is  Case  5.  This  data  was  collected 
and  analyzed  in  an  identical  fashion  as  the  previous  case,  but  instead  of  having  the 
simple  sinusoidal  pattern  as  seen  in  Case  3,  the  variance  in  the  geopotential  is  seen 
through  different  fluctuations  in  the  value  of  the  Z-component  of  the  inertial  angular 
momentum  but  the  constancy  of  this  component  remains  at  approximately  ICC5. 
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Specific  Angular  Momentum  -  J2,2 
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Figure  4.4:  Scenario  1,  Case  3:  Hk  versus  time 

Although  there  are  more  random-looking  variations  in  this  graph,  the  graph  continues 
to  be  cyclic.  However,  the  period  of  these  variations  appears  to  be  two-times  the 
period  of  the  J|  case.  This  implies  that  the  repetition  seen  in  the  previous  graph 
shows  the  variations  run  coincidentally  with  the  half-period  of  the  object’s  orbit  and 
the  repetition  in  Figure  4.5  is  coincident  with  the  period  of  the  orbit. 

The  last  case  analyzed  in  Scenario  1  is  the  Jl\  case.  The  constancy  of  the 
Z-component  of  the  inertial  angular  momentum  does  not  get  any  worse  with  the  in¬ 
creased  order  of  complexity  in  the  geopotential.  The  higher  order  terms  utilized  in 
Cases  3,  5,  7  and  9  limit  the  constancy  to  approximately  order  10~5  as  is  shown  in 
Figures  4.4,  4.5  and  4.6.  Additionally,  the  variations  seen  in  Figure  4.6  run  coinciden¬ 
tally  with  the  variations  in  Figure  4.5  validating  the  assumption  that  the  variations 
are  coincident  with  the  period  of  the  orbit.  Because  including  the  tesseral  and  sectoral 
harmonics  represent  a  more  realistic  geopotential  for  the  Earth’s  gravitational  field 
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Specific  Angular  Momentum  -  J4,4 
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Figure  4.5:  Scenario  1,  Case  5:  Hk  versus  time 

than  the  spherical  model  or  the  symmetric  model  developed  with  the  zonal  harmonics 
only,  10-5  is  the  the  level  of  constancy  that  can  be  expected  from  the  Z-component  of 
inertial  angular  momentum.  Graphs  from  Scenario  2  and  the  cases  not  cited  in  this 
section  can  be  found  in  Appendix  A. 

In  addition  to  graphing  the  results  from  the  Matlab®  program,  the  data  was 
imported  into  Excel®  and  analyzed  to  find  the  maximum  change  in  the  Z-component 
of  the  inertial  angular  momentum.  In  order  to  do  this,  within  Excel®,  a  column 
was  created  to  calculate  the  absolute  value  of  the  difference  between  Hf~{t)  —  Hk(t 0). 
Then  utilizing  the  built-in  function  in  Excel®  maximum  value  in  that  column  was 
determined.  The  results  for  all  four  scenarios  and  the  nine  cases  in  each  scenario  are 
summarized  in  Table  4.1. 
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Figure  4.6:  Scenario  1,  Case  9: 


Hk  versus  time 


Scenario  1 

Scenario  2 

Scenario  3 

Scenario  4 

Case  1 

2.74E-14 

1.23E-13 

3.88E-13* 

3.01E-13 

Case  2 

2.81E-13 

4.16E-13 

1.65E-12 

9.21E-13 

Case  3 

1.91E-05 

2.81E-05 

1.91E-05 

2.81E-05 

Case  4 

4.64E-13 

3.11E-13 

3.91E-12 

7.33E-13 

Case  5 

4.31E-05 

2.84E-05 

4.37E-05 

2.85E-05 

Case  6 

3.74E-13 

1.34E-13 

1.12E-12 

6.60  E-13 

Case  7 

4.50E-05 

2.85E-05 

4.63E-05 

2.87E-05 

Case  8 

2.90E-13 

1.13E-13 

1.31E-12 

2.63E-13 

Case  9 

4.48E-05 

2.97E-05 

4.54E-05 

3.00E-05 

Table  4.1:  Inertial  Angular  Momentum  Constancy  Results 
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V.  Right  Ascension  of  the  Ascending  Node 


5. 1  Introduction 

The  third  element  analyzed  for  constancy  is  the  time  rate  of  change  of  the 
right  ascension  of  the  ascending  node  (RAAN)  (fi).  Much  research  has  been  done 
on  the  effect  of  the  Earth’s  oblateness  term  in  the  geopotential,  J and  information 
can  be  found  in  several  different  books  and  papers.  The  author  chose  to  use  Dr. 
Wiesel’s  Modern  Astrodynamics  book  as  a  reference  on  this  topic.  In  the  two-body 
problem,  five  of  six  COEs  remain  constant  regardless  of  time,  but  as  stated  before,  the 
introduction  of  the  gravity  geopotential  limits  the  constancy  of  these  elements.  The 
effect  on  the  right  ascension  of  the  ascending  node  (0)  is  the  focus  of  this  chapter. 
In  his  book,  Dr.  Wiesel  addresses  the  rate  at  which  the  Node  changes  and  gives  an 
equation  in  the  “Lagrange  Planetary  Equations  Disturbing  Function  Form”  [10,  98]: 


dtt  _  1  d  R 

dt  na?  a/1  —  e2  sin  i  di 


(5.1) 


where  n  is  the  mean  motion  of  the  satellite,  a  is  the  semi-major  axis,  e  is  the  eccen¬ 
tricity  of  the  orbit  and  i  is  the  inclination. 

The  equation  for  D  is  not  useful  yet,  so  it  is  essential  to  determine  the  contri¬ 
butions  that  the  term  of  the  geopotential  adds  to  this  equation.  Therefore,  taking 
Equation  3.7  and  letting  n  —  2  and  rri  —  0  results  in  the  following  for  the  disturbing 
function: 

*3  =  -^ ^(WS-1)  (5.2) 

In  order  to  take  the  partial  of  the  disturbing  function  the  disturbing  function 
through  needs  to  be  written  in  terms  of  the  COEs  [10,  137].  Dr.  Wiesel’s  book 
contains  the  detailed  derivation  from  the  Lagrange  Planetary  Equation  form  to  the 
COE  form  [10,  137-140].  Rewriting  Equation  5.2  in  terms  of  the  classical  orbital 
elements  through  order  e2,  it  can  be  shown  that  [10,  140]: 

R ;  ■  0 11 + ?2)  (53) 
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Now  taking  the  partial  of  the  disturbing  function  in  Equation  5.3  with  respect  to 
inclination  as  detailed  in  Equation  5.1  the  following  is  an  expression  for  the  time  rate 
of  change  of  the  node  [10,  141]: 


O 


dCl  3nJ2Rl 

— —  = - 7j  cos  i 

dt  2 a°-  (1  -  e2)2 


(5.4) 


The  results  of  deriving  Equation  5.4  is  that  there  exists  a  real  value  that  can  be 
calculated  describing  how  changes  with  respect  to  time.  The  result  is  based  on 
the  term  in  the  geopotential  and  related  to  the  other  classical  orbital  elements. 
Although  other  equations  exist  applying  the  Brouwer-Lyddane  method  and  incorpo¬ 
rating  J2,  J3,  J4,  J5  and  ( J2 )2  terms,  Equation  5.4  is  sufficient  for  this  analysis  [3]. 
Additionally,  it  is  worth  stating  that  the  node  will  regress  for  prograde  orbits,  precess 
for  retrograde  orbits  and  is  unchanging  for  the  Molniya  orbit.  For  this  paper,  the 
example  scenarios  are  prograde  orbits  and  the  discussion  pertains  to  prograde  orbits. 


5.2  Analyzing  the  Constancy  of  the  Change  in  the  RAAN 

In  order  to  determine  the  constancy  in  the  time  rate  of  change  of  the  right 
ascension  of  the  ascending  node,  the  results  created  by  Dr.  Wiesel’s  obit  propagator 
serve  as  inputs  for  another  Matlab®  program  (See  Appendix  B).  The  same  four 
scenarios  were  tested  but  Case  1  was  not  analyzed  because  in  the  two-body  case  D 
does  not  change,  so  the  results  from  that  case  are  meaningless.  The  remaining  eight 
cases  (Case  2  -  Case  9)  as  stated  in  previous  chapters  were  analyzed  using  the  following 
initial  conditions: 

The  initial  conditions  for  Scenarios  1  and  3  are: 


Rijk  — 


t  1.05  \ 

0.0 

V  0-0  / 


{DU} 


(5.5) 
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r  0.0 

VIJK  =  0.845154254 

v  0.453464763 

The  initial  conditions  for  Scenarios  2  and  4  are: 

(  0.0 

Rijk  =  |  -1.15678606465 
0.0 

0.65744356375 

0.0 

0.65744356375 

The  constancy  of  the  time  rate  of  change  of  the  right  ascension  of  the  ascending 
node  is  evaluated  by  inputting  the  inertial  position  and  velocity  vectors  calculated  by 
Dr.  Wiesel’s  orbit  propagator  into  the  Matlab®  program  which  calculates  the  value 
of  the  right  ascension  of  the  ascending  node  at  each  time  step.  This  is  done  by  first 
calculating  the  Inertial  Angular  Momentum  using  Equation  4.1  and  then  calculating 
the  line  of  nodes  (n)  using  the  following  equation  [8,  118]: 

n  —  K  x  HjjK  (5-9) 

Using  the  fact  that  the  line  of  nodes  lies  in  the  Earth’s  equatorial  plane  we  know  the 
following  to  be  true  [9,  62]: 

77/  /s  /s 

n  =  7—^7  =  cosfll  +  sinflJ  (5.10) 

||n|| 

Applying  Equation  5.10  it  is  possible  to  solve  for  D  using  the  following  equation: 

n  -i  rij 

=  COS  7-^7 

|n| 
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(5.11) 


The  final  step  in  the  calculation  of  f2  is  to  do  a  quadrant  check  and  then  verify  that 
the  data  falls  in  a  good  range,  where  0  <  <  2n. 

The  Matlab®  program  calculates  for  each  time  step  and  outputs  the  data 
in  a  format  easily  imported  into  Excel®.  Once  the  data  is  imported  into  Excel®,  a 
graph  summarizing  the  data  was  created.  The  first  case  analyzed  is  Scenario  1  Case 
2,  which  includes  the  geopotential  expansion  through  the  term. 


Nodal  Regression  -  J2,0 


Figure  5.1:  Scenario  1,  Case  2:  Q  versus  time 

Figure  5.1  shows  that  the  time  rate  of  change  of  the  RAAN  is  very  close  to  being 
linear  with  time.  To  emphasize  this  point,  a  best-fit  trendline  is  plot  over  the  data. 
The  trendline  is  the  black  dotted  line  spanning  the  plotted  data.  The  equation  and 
the  value  for  R 2,  defined  as  “the  fraction  of  the  total  squared  error  that  is  explained 
by  the  model”,  are  also  noted  on  the  plot  [4],  The  R 2  value  is  close  to  1,  showing  how 
close  to  linear  is  with  respect  to  time  and  giving  a  sense  of  the  level  of  predictability 
for  this  data.  Additionally,  the  slope  of  the  line  is  the  value  determined  by  Equation 
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5.4.  In  order  to  determine  if  using  the  change  in  the  right  ascension  of  the  ascending 
node  produced  by  the  geopotential  expansion  through  the  term  is  appropriate,  the 
slope  from  Case  2  will  be  compared  with  the  slopes  from  the  remaining  cases  in  this 
scenario. 

The  next  case  analyzed  is  Case  3,  which  includes  zonal,  tessera!  and  sectoral 
harmonics  in  the  geopotential  expansion  through  J|.  This  case  allows  us  to  see  the 
differences  brought  about  by  including  the  tesseral  and  sectoral  harmonics  in  the 
geopotential  expansion.  Figure  5.2  shows  a  nearly  identical  regression  of  the  node 


Nodal  Regression  -  J2,2 


Figure  5.2:  Scenario  1,  Case  3:  versus  time 

with  respect  to  time.  The  slope  of  the  equation  defining  the  time  rate  of  change  of 
RAAN  is  constant  to  approximately  order  ICC8,  but  because  the  value  is  of  order  10”3 
the  result  is  constancy  to  five  significant  figures.  Adding  a  few  tesseral  and  sectoral 
harmonics  to  the  geopotential  did  not  seem  to  affect  the  linearity  of  as  shown  in 
the  R2  term.  From  the  change  in  the  slope,  the  effect  these  first  tesseral  and  sectoral 


30 


harmonic  terms  have  on  the  time  rate  of  change  of  RAAN  is  approximately  order 
ICC6.  This  implies  that  through  the  J|  term  in  the  expansion,  the  J3  derived  time 
rate  of  change  of  RAAN  is  an  adequate  estimation.  It  is  also  interesting  to  note  that 
the  R2  term  actually  increased  slightly  showing  that  this  case  is  more  linear  than  Case 
2. 

Noting  the  minimal  difference  in  the  slope  from  Case  2  to  Case  3,  analysis 
continued  with  the  next  case  evaluated  being  Case  4.  Case  4  includes  terms  through 
J4  in  the  geopotential  expansion.  This  case  gives  an  idea  of  how  well  the  J3  case  does 
with  the  inclusion  of  more  zonal,  tesseral  and  sectoral  harmonics  in  the  geopotential 
expansion.  Figure  5.3  shows  a  familiar  graph  to  the  first  two  but  this  time  the 


Nodal  Regression  -  J4,0 


Figure  5.3:  Scenario  1,  Case  4:  0  versus  time 

precision  of  the  slope  in  the  best-fit  equation  reduces  to  approximately  order  10-6. 
As  discussed  above,  because  the  original  value  is  of  order  10-3  the  result  is  accuracy  to 
three  significant  figures.  Once  again,  it  is  interesting  to  note  that  the  R 2  value  slightly 
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increased  getting  even  closer  to  1.  So  although  the  accuracy  of  the  assumption  is 
decreased  to  three  significant  figures,  the  results  are  showing  a  more  linear  regression 
in  the  node. 

The  final  case  analyzed  is  Case  9.  This  includes  the  full  expansion  of  the  geopo¬ 
tential  through  the  J2l  term.  Despite  the  significant  increase  in  terms  added,  the 
slope  of  the  best-fit  line  stayed  at  approximately  order  10-6.  This  result  shows  that 
there  is  no  impact  on  the  plausibility  of  using  the  assumption.  The  significant 
variance  came  from  adding  terms  through  Jj.  Any  further  expansion  of  the  geopo¬ 
tential  has  little  to  no  effect  on  the  accuracy  of  this  assumption.  Additionally,  the  R2 
value  increased  slightly  showing  that  for  this  orbit  and  cases  evaluated,  Case  9  is  the 
most  linear  case. 


Nodal  Regression  -  J21.21 


Figure  5.4:  Scenario  1,  Case  9:  Q  versus  time 

Although  there  is  not  a  full  description  of  all  the  scenarios  and  cases  in  this 
chapter,  a  summary  of  the  information  gathered  for  each  case  was  created.  In  order 


32 


to  analyze  the  data,  all  the  data  was  imported  into  Excel®  and  a  plot  was  created. 
On  each  one  of  the  graphs  there  is  a  best-fit  trendline  and  the  R 2  value  annotated. 
Tables  5.1  and  5.2  show  a  summary  of  the  results  for  all  relevant  cases  in  Scenarios 
1  and  2  respectively.  Additionally,  figures  are  available  in  Appendix  A  for  Scenario  2 
and  the  cases  not  specifically  addressed  in  this  chapter. 


Scenario  1 

Best-Fit  Trendline 

R2 

Case  2 

-0.00136047040768t  +  6.27690646984026 

0.984281416881 

Case  3 

-0.00136052157367t  +  6.27692733240726 

0.984282476559 

Case  4 

-0.00136344271778t  +  6.27692127388479 

0.984348728799 

Case  5 

-0.00136349596607t  +  6.27695859934003 

0.984349746613 

Case  6 

-0.00136394677915t  +  6.27691182676157 

0.984360161509 

Case  7 

-0.00136398333024t  +  6.27692569650047 

0.984360916666 

Case  8 

-0.00136408012100t  +  6.27690812910590 

0.984363189205 

Case  9 

-0.00136411270942t  +  6.27692461306990 

0.984363841352 

Table  5.1:  Scenario  1,  0  Results 


Scenario  2 

Best-Fit  Trendline 

R2 

Case  2 

-0.00069216082930t  +  4.71238962787044 

0.999999421426 

Case  3 

-0.00069213563793t  +  4.71237847388094 

0.999999420218 

Case  4 

-0.00069240615521t  +  4.71239121808520 

0.999999422356 

Case  5 

-0.00069237544806t  +  4.71235911459945 

0.999999419109 

Case  6 

-0.00069220307448t  +  4.71239091720158 

0.999999422679 

Case  7 

-0.00069216655269t  +  4.71236023876455 

0.999999419493 

Case  8 

-0.00069215279270t  +  4.71239087276285 

0.999999422789 

Case  9 

-0.00069211664459t  +  4.71235730338259 

0.999999419370 

Table  5.2:  Scenario  2,  0  Results 
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VI.  Data  Accuracy 


6. 1  Introduction 

Now  that  a  constancy  level  for  H,  H and  has  been  evaluated,  it  is  necessary 
to  investigate  the  levels  of  accuracy  for  these  constants  attainable  using  different 
sensors.  There  are  a  few  different  types  of  sensors  being  used  currently  to  track 
satellites  and  other  objects  orbiting  the  Earth.  The  most  popular  in  this  application 
are  optical,  radar  and  laser  instruments.  Different  types  of  sensors  have  varying 
limitations  based  on  the  technology  and  quality  of  equipment  installed  in  the  sensor. 
Vallado  provides  a  table  with  limitations  and  biases  of  United  States  owned  sensors, 
but  for  simplicity  the  author  will  evaluate  generalized  cases  with  no  ties  to  current 
sensors  to  get  an  idea  of  how  different  parameters  affect  the  quality  of  the  acquired 
data  [8,  245]. 

In  order  to  investigate  the  levels  of  accuracy  obtainable  by  different  sensors, 
position  and  velocity  data  will  be  simulated.  The  data  will  be  simulated  using  a 
simple  dynamics  model  and  the  amount  of  error  introduced  into  the  data  based  on 
observation  quality  and  quantity  will  be  determined.  Then  the  simulated  data  and 
observed  data  will  be  propagated  using  the  full  geopotential  through  the  Jfj1  term  to 
determine  the  correlation  with  observation  accuracy  and  the  accuracy  attainable  for 
the  Hamiltonian  Function,  Z-component  of  inertial  angular  momentum  and  the  time 
rate  of  change  of  the  right  ascension  of  the  ascending  node.  Finally,  the  constancy  in 
7i ,  Hk  and  D  attainable  utilizing  the  simulated  data  will  be  compared  with  previously 
determined  constancies  for  H,  H \  and  to  establish  the  basis  for  a  supplemental 
model  that  can  be  used  as  a  tool  for  data  mining. 

6.2  Simulating  R  and  V  Data 

In  order  to  determine  the  quality  of  data  produced  by  a  specific  sensor,  a  sim¬ 
plified  least  squares  estimator  was  written  in  Matlab®.  The  first  thing  required  in 
the  program  is  to  produce  data  to  observe.  Since  this  is  a  simplified  estimator,  a  flat- 
Earth  model  was  selected  to  represent  the  data.  The  data  was  created  using  basic 
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physics  kinematic  equations,  assuming  a  constant  acceleration: 


R(t)  —  Rq  +  Vot  +  -gt2  (6-1) 

V(t)  =  V0  +  gt  (6.2) 

where  R  is  the  position  vector  and  V  is  the  velocity  vector  at  a  given  time.  The 
subscript  “0”  indicates  the  initial  value  where  t  —  0,  and  g  is  the  acceleration  due 
to  gravity.  For  this  model,  the  assumed  value  of  g  is  the  local  gravity  at  the  surface 
of  the  Earth  or  9.81  ^  toward  the  center  of  the  Earth.  Using  Equation  6.1,  and 
inputting  the  following  initial  conditions: 


(6.3) 


(  0 

V0  =  6538.656905 

v  3600.9857675 

Simulated  data  was  produced  for  a  ten  minute  interval.  A  time  interval  of  ten  minutes 
was  selected  because  that  is  an  approximate  in- view  time  for  a  satellite  being  tracked 
by  a  scanning  sensor  on  Earth.  The  time  step,  At,  is  varied  from  0.01  s  to  60  s  and  is 
used  to  identify  the  impact  that  the  time  step  has  on  obtaining  a  good  estimate,  which 
will  be  addressed  later.  Now  in  order  to  make  this  data  look  more  like  actual  data, 
noise  was  incorporated  into  the  data.  Matlab®  has  a  built  in  function  that  creates 
random  numbers.  The  random  number  was  then  multiplied  by  the  value  sigma,  which 
represents  the  order  of  magnitude  of  the  noise  in  the  data,  and  subsequently  added  to 
the  simulated  data.  This  creates  more  realistic  data  which  is  then  incorporated  into 
least  squares  as  z,i,  the  observed  data. 
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Although  the  dynamics  model  is  not  linear,  using  the  observed  data,  it  is  possible 
to  run  linear  least  squares  to  get  an  approximate  initial  state  as  input  into  the  non¬ 
linear  least  squares  problem.  This  step  is  not  mandatory,  an  approximate  initial 
state,  or  an  educated  guess,  would  be  sufficient  input  into  the  non-linear  least  squares 
problem,  but  the  author  chose  to  get  an  estimate  through  the  linear  least  squares 
method.  Linear  least  squares  attempts  to  estimate  “the  state  of  a  linear  dynamical 
system  at  epoch  time”  [11,  67].  The  initial  state  is  given  by: 


X0  = 


(6.5) 


There  is  a  detailed  process  by  which  the  following  equations  are  derived  in  Dr.  Wiesel’s 
Modern  Orbit  Determination  book  but  these  derivations  will  not  be  discussed  in  this 
paper  [11,  67-70].  The  first  step  in  solving  for  the  initial  state  using  linear  least  squares 
is  to  find  the  “state  transition  matrix”  ($)  [11,  32],  The  state  transition  matrix  is  the 
term  that  transitions  the  state  from  one  time  to  a  different  time  and  can  be  calculated 
using  the  following  equation: 


$(Mo) 


dX(t) 

dX{t0) 


(6.6) 


Applying  Equations  6.1  and  6.2  the  state  as  a  function  of  time  is  given  by: 


X(t) 


Xo  +  x0t 
Vo  +  Vo  t 
zq  +  %o  t  +  \  g  t2 
x0 


yo 

\  Z0+gt  y 


(6.7) 
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Taking  the  partial  derivatives  in  Equation  6.6  gives  a  state  transition  matrix  as  follows: 


(t,  t0) 


(  1  0  0  t  0  0  ^ 
010010 
0  0  1  0  0  1 
0  0  0  1  0  0 
0  0  0  0  1  0 
v  0  0  0  0  0  1  } 


(6.8) 


The  next  step  is  to  define  a  matrix,  the  “observation  relation” ,  G(t)  that  predicts 
the  data  Zd(t)  [11,  75].  The  assumption  is  that  observation  data  is  the  inertial  position 
vector.  Therefore,  G(t)  is  given  by  the  following  equation: 


G{t) 


*  x(t)  ^ 

y(t) 

v  z(t)  / 


(6.9) 


Additionally,  the  H  matrix  is  the  linearization  of  the  observation  relation  and  can  be 
found  in  the  following  manner  [11,  76]: 


d  G(t) 
dX{t) 


(6.10) 


taking  the  partial  derivative  results  in  the  following  for  the  H  matrix: 


H  = 


^  1  0  0  0  0  0  ^ 
I  0  1  0  0  0  0 

v  0  0  1  0  0  0  y 


(6.11) 


Now,  for  short-hand,  we  will  define  the  T  matrix,  also  referred  to  as  the  “observation 
matrix”  [11,  67]: 

T(t)  =  H  $(Mo)  (6.12) 
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Then,  it  is  possible  to  define  the  “instrumental  covariance  matrix”  Q  as  equal  to 
sigma2,  where  sigma  is  the  order  of  magnitude  of  the  noise  in  the  data  due  to  the 
quality  of  the  instrument  [11,  68].  At  this  point,  applying  the  previous  equations 
makes  it  possible  to  solve  for  the  “state  covariance”  (Px)  defined  by  the  following 
equation  [11,  70]: 

Px  =  (6,13) 

The  covariance  is  calculated  by  first  running  a  for  loop  for  the  summation  and  then 
taking  the  inverse.  The  final  step  in  linear  least  squares  is  to  calculate  the  initial  state 
and  it  is  determined  by  the  following  equation  [11,  70]: 

N 

x0  =  pxJ2T?Qilz*  (6-14) 

i=  1 

Once  again  making  sure  that  the  summation  is  accomplished  in  the  afore  mentioned 
for  loop. 

Now  an  initial  guess  for  the  state  at  epoch  time  has  been  obtained  and  will  be 
applied  as  an  input  to  the  non-linear  least  squares  estimator.  Although  making  an 
educated  guess  would  make  it  possible  to  skip  the  above  process,  the  afore  mentioned 
process  provides  a  methodical  technique  for  obtaining  an  estimate  and  introduces 
many  of  the  terms  used  in  non-linear  least  squares.  Now  with  this  estimate  for  the 
initial  state,  an  update  of  the  state  is  required  and  can  be  found  using  the  following 
equation: 

X(t)  =  $(£,  to)  Xref(to)  +  Xp  (6.15) 

where  Xref(t0 )  is  the  estimate  calculated  in  the  previous  step  with  the  notation  up¬ 
dated  to  avoid  future  confusion  and  the  state  transition  matrix  is  the  same  as  defined 
in  Equation  6.8.  The  Xp  is  a  particular  solution  that  is  added  to  the  state  to  account 
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for  the  non-linear  terms  and  it  is  defined  below: 


0  1 

0 

\at2 

o 

o 

V  9t 


(6.16) 


Since  the  non-linear  terms  are  taken  care  of  by  adding  a  particular  solution  every 
time  the  state  is  updated,  the  remainder  of  the  problem  can  be  treated  as  a  linear 
problem  [11,  82-83].  The  next  step  is  to  calculate  observation  relation,  the  G  matrix, 
as  defined  in  Equation  6.9.  The  easiest  way  to  do  this  is  to  multiply  the  H  matrix 
by  the  updated  state  which  leaves  the  estimated  position  vector.  Then  calculating 
the  “residual  vector”  (f)  is  accomplished  using  the  relationship  between  the  observed 
data  and  the  observation  relation  matrix  [11,  77]: 


r(t)  =  zd(t )  -  G(t) 


(6.17) 


Now  using  Equation  6.12  the  T  matrix  is  calculated  and  it  is  used  to  calculate  the 
change  in  the  covariance  matrix  P$x  as  defined  in  Equation  6.18  [11,  77]: 

Psx  =  (TTQ~1Ty1  (6.18) 

Then  the  change  required  to  be  made  to  the  initial  state  estimate  is  calculated  to 
account  for  the  non-linear  terms  over-looked  in  the  original  estimate: 

SX0  =  P8xTTQ-]r  (6.19) 
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A  new  initial  state  is  then  calculated  using  the  following  equation: 

Xq  =  Xref(t0)  +  SX0  (6.20) 

Now  because  the  covariance  and  the  change  to  the  initial  state  needs  to  be  calculated 
for  every  At,  a  summation  is  calculated,  similar  to  what  was  done  in  the  linear  case, 
to  determine  the  total  values  of  P$x  and  TTQ~l,r  as  defined  in  Equations  6.18  and 
6.19.  Finally,  because  it  is  possible  to  gain  access  to  the  actual  initial  state,  the  error 
can  be  calculated  by  subtracting  the  initial  state  determined  through  least  squares 
from  the  initial  state  used  in  creating  the  simulated  data  (Equations  6.3  and  6.4).  The 
absolute  value  of  the  difference  between  these  terms  is  output  by  Matlab®  for  each  of 
the  position  and  velocity  components  as  well  as  the  magnitude  of  these  vectors  (see 
code  in  Appendix  B).  In  order  to  validate  the  computer  program,  the  noise  was  not 
included  with  the  data  and  executing  the  Matlab®  program  the  following  results  for 
the  error  were  obtained: 

Position  error  is:  0.000000000000e+000 
x-direction  0.000000000000e+000 
y-direction  4.758828049570e-012 
z-direction  6.984919309616e-010 
Velocity  error  is:  0.000000000000e+000 
x-direction  7.185018701990e-015 
y-direction  0.000000000000e+000 
z-direction  4.547473508865e-013 

Just  for  note,  the  program  was  run  with  At  equal  to  1  s  and  sigma  equal  to  1  m.  The 
results  from  this  run  of  the  program  show  that  the  error  between  the  calculated  initial 
state  and  the  initial  state  input  into  the  program  to  simulate  the  data  is  zero.  This 
validates  that  the  program  is  working  as  expected.  Now  that  a  computer  program 
has  been  written  to  simulate  the  data,  analysis  will  be  done  to  identify  the  correlation 
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between  observation  quality  and  quantity  and  the  error  in  the  position  and  velocity 
vectors. 

6.3  Quantifying  the  Data 

Utilizing  the  computer  code  for  the  least  squares  estimator  written  in  Matlab®, 
there  are  two  different  input  values  that  were  altered  to  assess  the  impact  on  the 
quality  of  the  estimate  for  the  initial  state  obtainable.  The  first  parameter  altered 
is  the  noise  level  in  the  data  and  it  is  based  on  the  quality  of  the  sensor  being  used, 
defined  by  the  instrumental  covariance  matrix.  In  order  to  test  the  impact  that  this 
parameter  has  on  the  quality  of  the  data,  values  from  0.1  m  to  1000  m  as  the  sigma  in 
increments  of  0.1  m  were  input  into  the  program.  The  second  parameter  altered  is  the 
time  step,  At.  For  each  increment  in  sigma  and  for  a  given  At,  the  program  outputs 
the  absolute  value  of  the  error  between  the  scalar  value  of  the  estimated  position  and 
velocity  vector  and  the  initial  position  and  velocity  vector  inserted  to  create  the  data. 
In  this  analysis,  time  steps  varying  from  0.05  s  up  to  60  s  were  implemented.  Below, 
Figure  6.1  shows  the  error  in  position  versus  sigma  for  the  shortest  time  step. 

In  this  analysis,  each  time  step  was  plotted  on  a  separate  graph  to  show  how 
position  error  changes  with  respect  to  the  accuracies  of  the  sensor  for  a  specific  time 
step.  The  first  case  where  the  time  step  is  very  small,  the  position  error  grows  linear 
with  respect  to  sigma.  A  trendline  is  plotted  over  the  data  and  the  best-fit  equation 
is  shown  on  the  graph  along  with  the  R2  term.  The  equation  shows  that  with  an  R2 
value  of  0.9995,  the  position  error  is  approximately  one-half  of  the  sigma  value.  This 
trend  continues  for  each  time  step  evaluated,  but  the  R 2  value  continually  decreases 
for  each  increase  in  the  time  step.  The  decrease  in  the  R2  value  implies  that  the 
data  is  more  scattered  providing  data  with  higher  variances  from  the  trendline.  For 
comparison,  the  position  error  versus  sigma  for  the  final  time  step,  where  At  is  equal 
to  60  seconds,  was  plotted.  Once  again,  Figure  6.2  shows  that  the  trend  for  the 
data  is  found  to  be  at  approximately  one-half  of  the  input  sigma  value.  However, 
because  the  R2  term  is  much  lower  than  it  was  in  Figure  6.1  there  are  more  data 
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Figure  6.1:  Initial  Position  Error  vs  Sigma  for  At  =  0.05s 

points  nearing  the  sigma  value.  Additionally,  an  analysis  was  done  to  calculate  the 
best-fit  trendline  and  the  R 2  value  for  each  time  step  and  the  results  for  all  the  cases 
are  summarized  in  Table  6.1.  The  graphs  for  time  steps  not  explicitly  discussed  in 
this  section  are  available  in  Appendix  A.  The  significance  of  this  analysis  is  that  for 


Time  (s) 

Best-Fit  Trendline 

R2 

0.05 

0.5001x  -  0.0321 

0.9995 

0.1 

0.4999x  +  0.0269 

0.9991 

0.5 

0.5003x  -  0.0911 

0.9957 

1 

0.4996x  +  0.1348 

0.9911 

5 

0.5003x  +  0.0737 

0.9571 

30 

0.4962x  +  0.6413 

0.8071 

60 

0.4984x  +  0.5105 

0.7078 

Table  6.1:  Initial  Position  Error  Results 


this  simplified  model  of  the  dynamics,  the  mean  value  of  the  data  is  on  the  order  of 
one-half  sigma.  This  means  that  for  a  sigma  value  of  one  meter,  the  error  in  position 
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Figure  6.2:  Initial  Position  Error  vs  Sigma  for  At  =  60s 

is  approximately  0.5  m.  Now,  because  position  is  on  the  order  of  106  m,  an  error 
of  0.5  m  provides  an  accuracy  in  position  to  seven  significant  figures  or  order  10~8. 
As  the  time  step  increases,  the  R 2  value  decreases  resulting  in  larger  variances  from 
the  best-fit  trendline.  However,  for  further  analysis,  the  values  as  they  appear  on  the 
best-fit  trendline  will  be  used  as  a  basis  for  the  accuracy  level.  Therefore,  the  range  of 
accuracies  in  the  position  vector  for  all  time  steps  goes  from  eight  significant  figures 
for  sigma  equal  to  0.1  m  down  to  four  significant  figures  for  sigma  equal  to  1000  m. 
Keep  in  mind  though,  that  as  the  time  step  increases  the  probability  that  the  error 
will  fall  on  or  near  the  best-fit  line  decreases,  as  evident  through  the  R 2  term.  At  a 
time  step  of  60  seconds  this  could  increase  the  error  to  the  order  of  the  sigma  value, 
which  would  give  a  range  of  accuracy  from  seven  to  three  significant  figures.  The 
accuracies  obtainable  using  the  mean  value  (value  obtained  from  the  trendline)  will 
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be  used  in  the  following  section  to  determine  the  impact  this  error  has  on  H *.  and 

a 

In  addition  to  position  error,  analysis  was  completed  to  determine  the  error  in 
the  velocity  vector.  Similarly  to  what  was  done  for  position,  the  absolute  value  of 
the  error  between  the  least  squares  calculated  velocity  and  the  original  velocity  vector 
used  to  create  the  simulated  data  was  calculated.  However,  a  slight  modification  to 
the  process  for  collecting  data  to  determine  the  mean  value  is  made  and  approximately 
100  data  points  for  every  magnitude  increase  in  sigma  will  be  evaluated.  Then,  as 
before,  the  time  step  is  varied  from  0.05  s  to  60  s  and  the  error  is  plotted  for  each 
time  step.  In  order  to  best  show  the  changes  in  error  as  sigma  increases,  the  data  is 
plotted  on  a  logarithmic  plot.  Figure  6.3  shows  the  first  time  step  analyzed. 


Figure  6.3:  Initial  Velocity  Error  vs  Sigma  for  At  =  0.05s 

Unfortunately,  the  plot  of  initial  velocity  error  versus  sigma  does  not  give  us  a 
strictly  linear  relation;  therefore,  a  power  law  equation  was  fit  to  the  data.  Although 
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the  K2  term  is  not  extremely  high,  the  power  law  fit  provides  a  nice  mean  value  to 
be  used  in  determining  the  accuracy  attainable  in  the  velocity.  Using  the  power  law 
fit  equation  plotted  in  Figure  6.3,  it  can  be  shown  that  the  velocity  error  ranges  from 
order  1CU'  —  at  sigma  equal  to  0.1  m  up  to  order  10-3  —  where  sigma  is  equal  to  1000 
m.  Because  the  velocity  is  on  the  order  of  103  — ,  accuracies  for  velocity  range  from 
ten  to  six  significant  figures.  Similar  to  the  position  error,  the  average  velocity  error 
increases  approximately  one  order  of  magnitude  for  every  order  of  magnitude  increase 
for  sigma.  However,  unlike  the  position  error  the  velocity  error  changes  with  respect 
to  the  At  term.  Using  a  time  step  of  two  times  the  first,  t  —  0.1  s,  results  in  an  order 
of  magnitude  increase  in  the  velocity  error  for  all  sigma  values.  Figure  6.4  shows  this 


Figure  6.4:  Initial  Velocity  Error  vs  Sigma  for  At  =  0.1s 

increase  in  the  velocity  error  and  the  power  law  fit  equation  can  be  used  to  calculate 
the  average  velocity  error  for  this  time  step.  The  velocity  error  starts  at  order  10-6 
—  for  a  sigma  value  of  0.1  m  and  increases  to  order  10~2  —  for  sigma  equal  to  1000 
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m.  The  basic  structure  of  the  plot  remains  the  same  showing  power  dependence  on 
sigma.  The  order  of  magnitude  error  in  velocity  seen  at  a  time  step  of  0.1  s  is  the 
same  for  At  values  of  0.5,  1,  and  5  seconds.  The  plots  for  each  of  these  cases  will  not 
be  shown  here  but  are  included  in  Appendix  A.  The  velocity  error  increases  again  by 
an  order  of  magnitude  for  the  30  and  60  second  time  steps. 
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Figure  6.5:  Initial  Velocity  Error  vs  Sigma  for  At  =  60s 

Figure  6.5  shows  the  initial  velocity  error  determined  using  least  squares  with 
a  time  step  of  60  seconds.  As  mentioned  previously,  the  error  for  this  time  step  is 
one  order  of  magnitude  higher  for  all  values  of  sigma  as  compared  with  At  equal  to 
0.1  s  and  two  orders  of  magnitude  higher  than  At  equal  to  0.05  s.  The  power  law  fit 
equation  for  this  data  was  used  to  calculate  the  range  of  velocity  error  expected  for 
sigma  values  from  0.1  to  1000  meters.  The  mean  value  of  the  initial  velocity  error 
is  calculated  to  be  of  order  10“5  —  at  the  lowest  sigma  value,  and  increases  to  order 
lO^1  —  at  the  highest  sigma  value.  Although  the  graph  for  the  30  second  time  step 
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is  not  included  in  this  section,  it  is  included  in  Appendix  A.  Additionally,  the  power 
law  fit  equation  and  R2  values  for  each  of  the  time  steps  analyzed  are  summarized  in 
Table  6.2. 


Time  (s) 

Power  Law  Fit  Equation 

R 2 

0.05 

0.00000793a:1  00506942 

0.85337492 

0.1 

0.00001 147a:0,99959821 

0.84789666 

0.5 

0.00002480a:1  00269564 

0.85068946 

1 

0.00003516a:1  00350272 

0.84703331 

5 

0.00007833a:1  00023009 

0.84223989 

30 

0. 00018685a:0"491486 

0.85111477 

60 

0. 00025541a:0'99164799 

0.84730417 

Table  6.2:  Initial  Velocity  Error  Results 


Now,  with  the  initial  position  and  velocity  errors  created  through  the  simplified 
least  squares  estimator,  it  is  possible  to  use  the  average  error  calculated  from  the 
best-fit  trendline  and  power  law  fit  equation  to  estimate  the  impact  these  errors  have 
on  7i,  Hk  and  Cl. 

6.4  Impact  on  the  Constants 

In  order  to  assess  the  impact  of  the  initial  position  and  velocity  error  on  Tt,  H k 
and  Cl,  the  initial  position  and  velocity  vectors  (Equations  6.3  and  6.4)  used  as  inputs 
to  create  the  data  were  converted  to  canonical  units  and  input  into  Dr.  Wiesel’s 
orbit  propagator  using  the  full  geopotential  expansion  through  the  J2\  term.  In  a 
similar  fashion  as  discussed  in  each  of  the  previous  chapters,  the  output  position  and 
velocity  vectors  from  the  propagator  were  used  to  calculate  the  expected  H,  Hk  and 
Cl  values  to  use  as  a  baseline.  Then  importing  the  data  into  Excel®,  baseline  values 
for  H,  Hk  and  Cl  were  included  in  separate  spreadsheets.  Next,  the  position  and 
velocity  vectors  calculated  using  the  Matlab®  least  squares  estimator,  with  errors 
approximately  equal  to  the  value  defined  by  the  associated  trendlines  (see  Tables  6.1 
and  6.2),  were  propagated  using  Dr.  Wiesel’s  orbit  propagator  incorporating  the  full 
geopotential  expansion  through  the  J2\  term.  Then  using  the  output  hies,  H,  Hk 
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and  £2  values  associated  with  these  initial  conditions  were  calculated  using  Matlab®. 
Since  the  author  is  interested  in  the  error  caused  in  these  three  constants  as  a  result 
of  the  error  in  the  initial  position  and  velocity  vectors,  the  new  values  for  and 

£7  were  input  into  the  respective  Excel®  spreadsheets  and  the  maximum  difference 
between  these  values  and  the  baseline  values  for  each  constant  over  the  ten  minute 
time  interval  was  calculated.  Analysis  for  all  time  steps  was  completed  to  include  the 
minimum  and  maximum  values  for  sigma  mentioned  previously.  Tables  6.3  through 
6.5  show  the  results  for  maximum  error  calculated  during  the  propagation  time  for 
the  Hamiltonian  Function,  Z-component  of  inertial  angular  momentum  and  the  time 
rate  of  change  of  the  right  ascension  of  the  ascending  node. 


0.1  m 

1000  m 

0.05  s 

6.76150E-09 

6.62014E-05 

0.1  s 

6.85740E-09 

6.81668E-05 

0.5  s 

6.85740E-09 

7.33783E-05 

1  s 

6.48999E-09 

6.87841E-05 

5  s 

6.93561E-09 

5.31149E-05 

30  s 

9.34998E-09 

5.54704E-05 

60  s 

1.16882E-08 

8.75611E-05 

Table  6.3:  Hamiltonian  Function  Impacts 


0.1  m 

1000  m 

0.05  s 

9.21223E-09 

8.07959E-05 

0.1  s 

9.04901E-09 

8.33492E-05 

0.5  s 

9.04901E-09 

8.59880E-05 

1  s 

8.22583E-09 

8.38054E-05 

5  s 

8.31630E-09 

6.38745E-05 

30  s 

9.50784E-09 

8.00142E-05 

60  s 

8.29201E-09 

1.40391E-04 

Table  6.4:  Z-Component  of  Inertial  Angular  Momentum  Impacts 


The  error  in  the  Hamiltonian  Function  and  the  Z-component  of  inertial  angular 
momentum  are  on  the  same  order  of  magnitude  as  the  error  in  the  initial  position 
vectors.  As  mentioned  earlier,  initial  position  error  was  on  the  order  of  one-half 
sigma  which  gives  us  a  range  from  eight  significant  figures  for  the  minimum  value 
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0.1  m 

1000  m 

0.05  s 

2.40150E-07 

2.30941E-03 

0.1  s 

2.44420E-07 

2.38862E-03 

0.5  s 

2.44420E-07 

2.58224E-03 

1  s 

2.31306E-07 

2.41162E-03 

5  s 

2.49130E-07 

1.85890E-03 

30  s 

3.44855E-07 

1.84797E-03 

60  s 

4.43634E-07 

3.03183E-03 

Table  6.5:  Nodal  Regression  Impacts 


of  sigma  analyzed  down  to  four  significant  figures  for  the  maximum  value  of  sigma. 
The  Hamiltonian  Function  and  Z-component  of  inertial  angular  momentum  follow 
the  same  pattern  showing  that  the  accuracy  of  these  constants  is  directly  related 
to  the  accuracy  of  the  position  attainable,  which  is  related  to  the  accuracy  of  the 
sensor.  The  time  step  does  not  appear  to  be  directly  related  to  the  errors  in  7 H.  and 
Hk  but  that  is  a  result  of  selecting  data  near  the  trendline.  As  the  time  step  increases, 
the  probability  that  these  values  will  be  on  or  near  the  trendline  decreases,  possibly 
impacting  the  accuracies  up  to  one  order  of  magnitude  as  stated  previously. 

The  error  in  the  time  rate  of  change  of  the  right  ascension  of  the  ascending  node 
appears  different  from  the  others  because  it  is  not  solely  dependent  on  the  position 
accuracy.  With  an  error  starting  at  order  10-'  for  the  minimum  value  of  sigma  and 
increasing  to  order  10-3  for  the  maximum  value  of  sigma  is  two  orders  of  magnitude 
lower  than  would  be  expected  if  the  error  was  limited  by  the  accuracy  of  the  position 
alone.  The  data  created  for  the  nodal  regression  shows  that  initially  the  error  is  on 
the  order  of  10~9  and  10“5  for  sigma  values  of  0.1  and  1000  m  respectively,  but  as 
the  propagation  time  progresses  the  error  increases  showing  the  dependence  on  time 
as  well.  Therefore,  analysis  was  done  to  find  the  slope  related  to  the  regression  of  the 
node  for  each  case.  Then  the  absolute  value  of  the  difference  between  the  slope  of  the 
baseline  case  with  each  of  the  remaining  cases  was  determined.  The  difference  in  the 
slope  for  each  case  analyzed  is  summarized  in  Table  6.6.  Therefore  the  combination  of 
the  initial  error  in  the  position  and  the  time  period  for  which  the  data  was  propagated 
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0.1  m 

1000  m 

0.05  s 

7.7701E-11 

7.57278E-07 

0.1  s 

7.8852E-11 

7.83315E-07 

0.5  s 

7.8852E-11 

8.47747E-07 

1  s 

7.465E-11 

7.92397E-07 

5  s 

8.1049E-11 

6.05814E-07 

30  s 

1.11668E-10 

6.11538E-07 

60  s 

1.43509E-10 

9.87744E-07 

Table  6.6:  Nodal  Regression  Slope  Impacts 

contributed  to  the  error  in  the  regression  of  the  RAAN  as  seen  in  Table  6.5.  With 
this  information  it  is  now  possible  to  look  at  what  level  of  fidelity  of  the  geopotential 
expansion  is  appropriate  for  application  in  a  supplemental  model. 

6.5  A  Supplemental  Model 

According  to  the  level  of  accuracies  attained  in  the  Hamiltonian  Function  by 
using  simulated  data  it  does  not  appear  necessary  to  use  geopotential  terms  through 
Jf*.  Comparing  the  order  of  magnitude  in  Table  6.3  with  Table  3.2,  geopotential 
expansion  through  J. ®  is  sufficient  for  sensors  with  sigma  on  the  order  of  1000  m.  This 
is  because  the  error  seen  in  the  Hamiltonian  is  on  the  order  of  10~5  which  matches 
the  level  of  constancy  found  when  evaluating  the  Hamiltonian  Function  through  J9. 
For  a  slightly  better  sensor  where  sigma  is  on  the  order  of  100  m,  an  increase  in  the 
geopotential  terms  through  J|  is  appropriate.  Table  3.2  shows  that  the  constancy 
attained  in  the  Hamiltonian  Function  with  geopotential  expansion  through  this  term 
was  approximately  order  10-6,  which  matches  the  error  level  for  a  100  m  sensor.  In 
order  to  get  accuracies  of  10-9  more  cases  similar  to  those  listed  in  Table  3.2  were 
analyzed  and  it  was  found  that  just  decreasing  the  geopotential  expansion  to  J|q 
results  in  accuracies  on  the  order  of  10~7  which  would  not  be  sufficient  to  observe 
order  10-9  errors  in  the  Hamiltonian  Function.  Therefore,  sensors  attaining  meter 
accuracy  or  better  require  a  model  incorporating  a  geopotential  expansion  through 
the  J21  term. 
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The  Z-component  of  inertial  angular  momentum  results  displayed  in  Table  6.4 
show  that  for  higher  sigma  values,  the  error  is  on  the  order  of  10-5  which  coincides 
with  the  highest  level  of  constancy  that  can  be  expected  for  H k.  Table  4.1  shows  that 
for  any  geopotential  (ignoring  the  zonal  harmonics  only  cases)  the  highest  level  of 
constancy  is  of  order  10-5  when  using  geopotential  terms  from  J|  all  the  way  through 
A\-  Therefore,  it  is  not  necessary  to  include  any  more  terms  in  the  geopotential 
expansion  beyond  the  J|  term. 

Finally,  looking  at  the  results  from  the  time  rate  of  change  of  the  right  ascension 
of  the  ascending  node,  the  change  in  the  slope  shown  in  Table  6.6  is  so  low,  the 
difference  would  not  be  noticeable  in  the  model  as  discussed  in  Chapter  V.  Table  5.1 
shows  that  the  maximum  difference  between  slopes  is  on  the  order  of  10~6.  Therefore, 
the  minor  differences  brought  about  by  this  particular  data  set  are  not  the  limiting 
factor  in  the  accuracy  attainable  in  the  time  rate  of  change  of  the  right  ascension  of 
the  ascending  node  and  the  model  should  be  developed  based  on  the  highest  level 
of  accuracy  attainable  in  the  slope  using  the  lowest  possible  geopotential  expansion. 
As  previously  stated,  the  largest  variance  occurs  between  A  and  and  adding 
geopotential  terms  beyond  this  term  results  in  no  additional  variance.  Therefore,  a 
model  utilizing  the  geopotential  expansion  through  the  A  is  appropriate.  Table  6.7 
summarizes  the  level  of  geopotential  expansion  appropriate  for  different  sigma  values. 
Current  systems  have  sigma  values  on  the  order  of  hundreds  of  meters,  implying  that 


Meter  or  Better 

Tens  to  Hundreds  of  Meters 

Thousands  of  Meters 

H 

j-21 

J21 

74 

J  4 

A 

Hk 

A 

A 

A 

0 

A 

A 

A 

Table  6.7:  Summary  of  Supplemental  Model 


the  geopotential  expansion  appropriate  with  this  order  of  sigma  would  be  sufficient 
for  implementation  in  data  mining  techniques. 
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VII.  Conclusion 


7.1  A  Growing  Problem 

As  the  number  of  objects  orbiting  the  Earth  continue  to  increase  through  cur¬ 
rent  launch  operations,  failing  satellites  and  a  decrease  in  the  size  of  objects  being 
tracked,  the  probability  that  an  object  will  be  input  into  the  Satellite  Catalog  multiple 
times  increases.  Additionally,  because  tracking  satellites  requires  high-level  precision, 
extraneous  data  in  the  SATCAT  is  highly  undesirable.  In  order  to  limit  the  number 
of  uncorrelated  tracks  being  added  to  the  SATCAT,  it  is  beneficial  to  supplement  the 
current  system  for  tracking  and  updating  orbits  with  a  model  that  uses  elements  of 
the  motion  that  remain  more  constant  than  the  classical  orbital  elements. 

7.2  Analysis  Summary 

In  order  to  develop  an  appropriate  supplemental  model,  research  was  conducted 
to  determine  the  constancy  of  three  elements  of  the  motion,  the  Hamiltonian  Func¬ 
tion,  Z-component  of  angular  momentum  and  the  time  rate  of  change  of  the  right 
ascension  of  the  ascending  node.  In  addition  to  the  investigation  of  the  constancy  at¬ 
tainable  by  these  three  elements,  simulated  data  was  analyzed  to  identify  accuracies 
attainable  through  current  estimation  procedures  based  on  the  quality  and  quantity 
of  the  observations. 

The  Hamiltonian  Function  was  shown  to  provide  accuracies  on  the  order  of 
twelve  significant  figures.  Additionally,  the  simulated  data  showed  that  using  a  sensor 
with  accuracies  of  meter  level  or  better  would  require  using  a  geopotential  expansion 
through  the  Jfj1  term.  Sensors  with  accuracies  on  the  order  of  ten  meters  to  hundreds 
of  meters  would  require  a  model  using  a  geopotential  expansion  through  the  Jf  term, 
while  a  sensor  on  the  order  of  1000  meters  would  only  require  expansion  through  the 
J 2  term.  Therefore,  depending  on  the  accuracy  of  the  sensor  in  use,  the  fidelity  of 
the  model  will  need  to  be  adjusted. 

The  Z-component  of  inertial  angular  momentum  is  another  element  that  was 
evaluated  to  determine  constancy.  Using  several  zonal  and  sectoral  cases,  the  con- 
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stancy  of  this  element  was  shown  to  be  on  the  order  of  four  significant  figures  with 
the  inclusion  of  the  tessera!  and  sectoral  harmonics.  The  theoretical  level  of  accu¬ 
racy  coincides  with  the  level  of  accuracy  attainable  using  a  sensor  with  accuracy  on 
the  order  of  1000  meters.  Therefore,  the  model  should  use  the  geopotential  expansion 
through  the  J|  term  in  order  to  attain  the  best  available  accuracy  in  the  Z-component 
of  inertial  angular  momentum. 

The  final  element  of  the  motion  evaluated  for  accuracy  is  the  time  rate  of  change 
of  the  right  ascension  of  the  ascending  node.  The  right  ascension  of  the  ascending 
node  decreases  linearly  with  time  at  a  rate  dependent  on  the  classical  orbital  elements. 
The  most  significant  variance  in  U  occurs  with  the  expansion  through  the  term 
resulting  in  accuracies  of  approximately  three  significant  figures.  The  simulated  data 
was  producing  errors  of  much  lower  values  than  the  model  at  this  geopotential  ex¬ 
pansion  will  show.  Therefore,  an  appropriate  model  is  one  utilizing  the  geopotential 
expansion  through  the  term. 

7.3  A  Model  for  Improvement 

Each  of  the  elements  of  the  motion  evaluated  for  constancy  require  a  slightly 
different  level  of  expansion  for  the  geopotential  to  be  used  in  the  model.  Utilization  of 
this  model  for  data  mining,  along  with  the  current  method  has  potential  to  enhance 
the  level  of  accuracy  of  objects  being  tracked  by  introducing  elements  of  the  motion 
that  remain  constant  for  longer  periods  of  time.  An  increase  in  tracking  accuracy 
will  decrease  the  number  of  uncorrelated  tracks  appearing  in  the  Satellite  Catalog, 
thereby  improving  the  available  data  and  enhancing  awareness  of  the  objects  orbiting 
the  Earth. 
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Hamiltonian  [DUA2mJA2] 


Appendix  A.  Additional  Figures 


A.l  Scenario  1:  Hamiltonian  Function  Graphs 

Hamiltonian  Function  -  J2,2 


Figure  A.l:  Scenario  1,  Case  3:  7i  versus  time 
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Hamiltonian  [DUA2/TUA2]  Hamiltonian  [DUA2mJA2] 


Hamiltonian  Function  -  J4,0 


-0.54510066632180 

-0.54510066632185 

-0.54510066632190 

-0.54510066632195 

-0.54510066632200 

-0.54510066632205 

-0.54510066632210 

-0.54510066632215 

-0.54510066632220 

-0.54510066632225 

-0.54510066632230 

-0.54510066632235 

Time  [TU] 

Figure  A. 2:  Scenario  1,  Case  4:  TL  versus  time 


Hamiltonian  Function  -  J4,4 


Figure  A. 3:  Scenario  1,  Case  5:  H  versus  time 
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Hamiltonian  [DUA2/TUA2]  Hamiltonian  [DUA2/TUA2] 


Hamiltonian  Function  -  J8,0 


Figure  A. 4:  Scenario  1,  Case  6:  7i  versus  time 


Hamiltonian  Function  -  J8,8 


Figure  A. 5:  Scenario  1,  Case  7:  H  versus  time 
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Iz.ru/?  *nal  ubkioimuibh 


Hamiltonian  Function  -  J21,0 


Figure  A. 6:  Scenario  1,  Case  8:  7i  versus  time 
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Z -Component  [DUA2/TU] 


A. 2  Scenario  1:  Angular  Momentum  Graphs 

Specific  Angular  Momentum  -  J4,0 


Figure  A. 7:  Scenario  1,  Case  4:  versus  time 
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Z -Component  [DUA2/TU]  Z -Component  [DUA2/TU] 


Specific  Angular  Momentum  -  J8,0 


Figure  A. 8:  Scenario  1,  Case  6:  Hk  versus  time 


Specific  Angular  Momentum  -  J8,8 


Figure  A. 9:  Scenario  1,  Case  7:  Hk  versus  time 
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Node  [rad] 


A. 3  Scenario  1:  Nodal  Regression  Graphs 


Nodal  Regression  -  J4,4 


Figure  A. 10: 


Scenario  1,  Case 


5:  versus  time 
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Node  [rad]  Node  [rad] 


Nodal  Regression  -  J8,0 


Figure  A. 11:  Scenario  1,  Case  6:  12  versus  time 


Nodal  Regression  -  J8,8 


Figure  A. 12:  Scenario  1,  Case  7:  versus  time 
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Node  [rad] 


Nodal  Regression  -  J21.0 


Figure  A. 13:  Scenario  1,  Case  8:  versus  time 
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Hamiltonian  [DUA2/TUA2] 


A. 4  Scenario  2:  Hamiltonian  Function  Graphs 


Hamiltonian  Function  -  J0,0 


Figure  A.  14:  Scenario  2,  Case  1:  7i  versus  time 
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Hamiltonian  [DUA2mJA2]  Hamiltonian  [DUA2mJA2] 


Hamiltonian  Function  -  J2,0 


Figure  A.  15:  Scenario  2,  Case  2:  Tt  versus  time 


Hamiltonian  Function  -  J2,2 


Figure  A.  16:  Scenario  2,  Case  3:  1H.  versus  time 
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Hamiltonian  [DUA2mJA2]  Hamiltonian  [DUA27TUA2] 


Hamiltonian  Function  -  J4,0 


Figure  A. 17:  Scenario  2,  Case  4:  Ti  versus  time 


Hamiltonian  Function  -  J4,4 


Figure  A.  18:  Scenario  2,  Case  5:  H  versus  time 


65 


Hamiltonian  [DUA27TUA2]  Hamiltonian  [DUA2mJA2] 


Hamiltonian  Function  -  J8,0 


Figure  A.  19:  Scenario  2,  Case  6:  7i  versus  time 


Hamiltonian  Function  -  J8,8 


Figure  A. 20:  Scenario  2,  Case  7:  H  versus  time 
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Hamiltonian  [DUA2mJA2]  Hamiltonian  [DUA27TUA2] 


Hamiltonian  Function  -  J21,0 


Figure  A. 21:  Scenario  2,  Case  8:  7i  versus  time 


Hamiltonian  Function  -  J21.21 


Figure  A. 22:  Scenario  2,  Case  9:  Tt  versus  time 
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Z -Component  [DUA2/TU] 


A. 5  Scenario  2:  Angular  Momentum  Graphs 


Specific  Angular  Momentum  -  J0,0 


Figure  A. 23:  Scenario  2,  Case  1:  Hk  versus  time 
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Z -Component  [DUA2/TU]  Z -Component  [DUA2/TU] 


Specific  Angular  Momentum  -  J2,0 


Figure  A. 24:  Scenario  2,  Case  2:  Hk  versus  time 


Specific  Angular  Momentum  -  J2,2 


Figure  A. 25:  Scenario  2,  Case  3:  Hk  versus  time 
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Z -Component  [DUA2/TU]  Z -Component  [DUA2/TU] 


Specific  Angular  Momentum  -  J4,0 


Figure  A. 26:  Scenario  2,  Case  4:  Hk  versus  time 


Specific  Angular  Momentum  -  J4,4 


Figure  A. 27:  Scenario  2,  Case  5:  Hk  versus  time 
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Z -Component  [DUA2/TU]  Z -Component  [DUA2/TU] 


Specific  Angular  Momentum  -  J8,0 


Figure  A. 28:  Scenario  2,  Case  6:  Hk  versus  time 


Specific  Angular  Momentum  -  J8,8 


Figure  A. 29:  Scenario  2,  Case  7:  Hk  versus  time 
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Z -Component  [DUA2/TU]  Z -Component  [DUA2/TU] 


Specific  Angular  Momentum  -  J21 ,0 


Figure  A. 30:  Scenario  2,  Case  8:  Hk  versus  time 


Specific  Angular  Momentum  -  J21,21 
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Figure  A. 31:  Scenario  2,  Case  9:  Hk  versus  time 
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Node  [rad] 


A. 6  Scenario  2:  Nodal  Regression  Graphs 


Nodal  Regression  -  J2,0 


Figure  A. 32:  Scenario  2,  Case  2:  versus  time 
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Node  [rad]  Node  [rad] 


Nodal  Regression  -  J2,2 


Figure  A. 33:  Scenario  2,  Case  3:  versus  time 


Nodal  Regression  -  J4,0 


Figure  A. 34:  Scenario  2,  Case  4:  versus  time 


74 


Node  [rad]  Node  [rad] 


Nodal  Regression  -  J4,4 


Figure  A. 35:  Scenario  2,  Case  5:  versus  time 


Nodal  Regression  -  J8,0 


Figure  A. 36:  Scenario  2,  Case  6:  versus  time 
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Node  [rad]  Node  [rad] 


Nodal  Regression  -  J8,8 


Figure  A. 37:  Scenario  2,  Case  7:  versus  time 


Nodal  Regression  -  J21,0 


Figure  A. 38:  Scenario  2,  Case  8:  versus  time 
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Node  [rad] 


Nodal  Regression  -  J21.21 


Figure  A. 39:  Scenario  2,  Case  9:  versus  time 
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A. 7  Quantifying  the  Data:  Initial  Position  and  Velocity  Error  Graphs 


Figure  A. 40:  Initial  Position  Error  vs  Sigma  for  At  =  0.1s 
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Position  Error  [m]  Position  Error  [m] 
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Sigma  [ni] 

Figure  A. 41:  Initial  Position  Error  vs  Sigma  for  At  =  0.5s 


Figure  A. 42:  Initial  Position  Error  vs  Sigma  for  At  =  Is 
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Position  Error  [m]  Position  Error  [m] 
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Figure  A. 43:  Initial  Position  Error  vs  Sigma  for  At  =  5 s 
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Figure  A. 44:  Initial  Position  Error  vs  Sigma  for  At  =  30s 
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Velocity  Error  [m/s]  Velocity  Error  [m/s] 
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Figure  A. 45:  Initial  Velocity  Error  vs  Sigma  for  At  =  0.5s 


Figure  A. 46:  Initial  Velocity  Error  vs  Sigma  for  At  =  Is 
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Velocity  Error  [m/s]  Velocity  Error  [m/s] 
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Figure  A. 47:  Initial  Velocity  Error  vs  Sigma  for  At  =  5s 
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Figure  A. 48:  Initial  Velocity  Error  vs  Sigma  for  At  =  30s 


82 


1 

6 

11 

16 

21 

26 

31 

36 

41 

46 


Appendix  B.  Matlab  Code 


B.l  Hamiltonian  Function  Matlab  Code 


Listing  B.l:  Hamiltonian  Function 

(app  endix2  /  Tot  alEnergy  auto .  m ) 

function  [TE] =TotalEnergy_auto (C , S) 

7. 7. 7. 7. 1 1 1 1 1 1 1 1 7. 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 7. 7. 

7.  This  function  takes  the  coefficients  C  and  S 
7.  created  from  the  EGM96  function  and  using  the 
7.  inertial  position  and  velocity  of  a  satellite 
°l  computes  the  Hamiltonian  of  the  system 

7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 


clc 

7.  Define  constants 

7.w_Earth  =  [0  ;  0  ;  0 . 000072921 1585530d0]  ; 

7.  From  Vallado’s  book 

w_Earth  =  [0 ; 0 ; 0 . 058833592d0] ; 

7.  Canonical  Units  From  Wiesel 

ord  =  l;  7«must  be  at  least  1 

7oOrder  1  gets  me  COO,  S00  and  POO 
maxtess  =0 ; 

7.  Get  the  data  from  the  file 

In=dlmread (’I:\My  Documents \Thesis\EarthSat\Data\.  .  . 

Circular45ShortTBP\ecr . txt ’ ) ; 
limit=length ( In) ; 

7«limit=30  ; 


7.  Put  the  data  into  new  matrices  time,  x  and  xdot 
a=l  ; 
b  =  1 ; 

while  b  <  limit 

time (a) =In (b) ; 
a=a  + 1  ; 
b=b+3 ; 

end 

c  =  l ; 
d  =  2  ; 

while  d  <  limit 

x (c , 1) =In (d  ,  1)  ; 
x (c , 2) =In (d  ,  2)  ; 
x  (c  ,  3)  =In (d  ,  3)  ; 
c  =  c  +  l  ; 
d=d+3 ; 

end 
e  =  l  ; 
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f=3; 

while  f  <  limit+1 

xdot(e,l)=In(f,l) 
xdot (e ,2) =In(f ,2) 
xdot (e ,3) =  In  (f ,3) 
e  =  e  + 1 ; 
f=f +3; 

end 


56  7,  Call  the  subroutine  EGM96  to  produce  C  and  S 
[C , S , mu , R_Earth]  =  EGM96(ord); 

C=double (C) ; 

S=double (S) ; 
mu=double  (1)  ; 

61  R_Earth=double ( 1 ) ; 

%  Input  the  position  and  velocity  vectors 
limit2=limit/3; 


66  /  Open  the  files  for  writing 

fid  =  fopen( 1  I : \ My  Document  s\Thesis\EarthSat\Data\Circular45ShortTBP 
\TotalEnergy_grob.txt1,  ’wt’); 
fid2  =  fopen( ’I :\My  Document s\Thesis\EarthSat\Data\.  .  . 

Circular45ShortTBP\TotalEnergy_grob_graph . txt ’ ,  ’ wt ’ ) ; 

for  y=l : limit2 

xl =  [x (y  ,  1 )  x (y  ,  2)  x(y,3)]’; 

71  xdotl=[xdot(y,l)  xdot(y,2)  xdot(y,3)]’; 

7.  Calculate  the  inertial  velocity 

crossterm =double (w_Earth* ( (xl (2 , 1) *  xdot 1 (1 , 1) ) -(xl (1 , 1) *xdotl . 

(2,1))))  ; 

crossterm  =  crossterm  (3)  ; 

KE1  =0 . 5*  dot  (  xdot  1  ,  xdot  1  )  +  cr  os  st  erm  ;  /  (m/sec)~2 


76 


81 


7.  Calculate  the  potential  energy 
rl=norm(xl);  7.  m 
r 1 =double ( r 1 ) ; 

U_Coef  f  1  =-mu/rl  ;  7.  (m/  sec)  '2 
U_Coeffl=double(U_Coeffl) ; 
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7.  Calculate  the  summat  ion  for  Geopotential 
U1 ( 1 ) =1 ; 

U1 =double (Ul) ; 

7.  Create  a  matrix  P  for  the  Legendre  Polynomials 
7.  where  Pnm  is  actually  P(m  +  l,n  +  l) 

PI (1 , 1) =1 ; 


91  7.  For  shorthand  let  cos  (theta)  =  p 

7,  In  ECI  coordinate  frame  cos  (theta) 

pl=(xl(3)/rl); 

pi =double ( pi ) ; 

7.f  printf  (  ’  pi  =  7.0 . 15f  ’  ,pl)  ; 


z/r 
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96 

101 

106 

111 

116 

121 

126 

131 

136 

141 


PI (1 ,2) =(pl)  ; 

PI (2 ,2)=-Pl (1 , l)/sqrt (l-pl-2)  ; 

'/.  for  m  =  0  see  equation  4.44  page  113  in  Modern  Astrodynamics 
for  n=2 : ord 

Pl(l,n+l)=(((2*(n-l)+l)*pl*Pl(l,n))-((n-l)*Pl(l,n-l)))/(n)... 

) 

f or  m  =  0 : n-2 

PI (m  +  2,n)  =  (((n-l) -m)*pl*Pl (m+1 ,n) -((n-l)+m)*Pl (m  +  1 ,n.  .  . 
-1) ) / sqrt (1-pl ~2) ; 

end 

end 

UE1 (l)=U_Coeffl*Ul (1) ; 

TE1 (1) =KE1+UE1 (1) ; 
i  =  2; 

for  n=l : ord -1 

if  maxtess  <  n 

mmax=maxtess ; 

else 

mmax=n ; 

end 

for  m=0 : mmax 

ang  =  at an 2 (xl(2,l),xl(l,l)); 

inter  =  C(n  +  l ,m  +  l) *  cos (m*  ang )+S(n  +  l ,m  +  l) *  sin (m*  ang )  ; 

U1 (i)=Ul ( i  —  1 ) + ( PI (m  +  1 , n  +  1) *  inter *(( rl /R_Earth ) ~ (-n) ) )  ; 
Ul=double (U1 ) ; 

UE1 (i)=U_Coef f 1*U1 (i) ; 

UEl=double (UE1 ) ; 

TE1 ( i ) =KE1 +UE1 (i) ; 

TEl=double ( TE1 ) ; 

TED IFF  1 (i)=TEl (i)-TEl ( i  —  1 )  ; 
i=i+l ; 
end 

end 

fprintf (fid  ,  ’Time :  %1 . 6f \n ’  ,time (y) )  ; 

k  =  l ; 

for  n=0 : ord - 1 

if  maxtess  <  n 

mmax=maxtess ; 

else 

mmax=n ; 

end 

for  m=0 : mmax 
if  n==ord-l 
if  m==n 

fprintf  (fid,  ’Order  J“/01  .  Of  ,  %1 . Of  :\n’  ,n,m)  ; 
fprintf (fid ,’ Total  Energy  is  %0 . 15f \n  ’  ,  TE1 (k) )  ; 

end 

end 

if  n==ord-l 
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146 


7«0 . 15f  \n  ’  ,time(y)  ,TEl(k 


151 


156 


161 


if  m==mmax 

fprintf (f id2  ,  ’ 7.4 . 14f 

))  ; 

end 
end 

7,  i  f  n  =  =  m 

7.  if  n  =  =ord-l 

•/,  fprintf  (fid,’ 7. 1.14f  "/oO.lSfXn’.time 

(y)  ,  TE1 (k)  )  ; 

7,  end 

7.  end 

k=k+ 1 ; 

end 

end 

fprintf (fid,  ’  \n ’ )  ; 


end 

fclose (fid) ; 
fclose (fid2) ; 
fprintf ( ’ D0NE\n  ’  ) ; 
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Listing  B.2:  EGM  96  Data 

(appendix2/EGM96.m) 


function  [C , S ,mu , R_Earth] =EGM96 (ord) 


7. 7. 1 7. 1 1 1 1 1 1 1 7. 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 7. 1 1 1 1 1 1 1 1 1 1 7. 7.  %  7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 


7.  This  function  takes  data  from  the  gravity  model  in 
7.  egm96.dat  and  creates  two  orderXorder  matrices  C  and  S 
7.  that  coordinate  with  the  normalized  coefficients  Cnm 
7.  and  Snm  .  The  matrix  is  tied  to  the  coefficient  in  the 
7.  following  way:  Cnm  =  C(n+l,m  +  l)  and  Snm  =  S  (n  +  1  ,  m  + 1 ) 


7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 


7.  Initialize  the  C  and  S  matrices 
C=zeros (ord , ord) ; 

S  =  zeros (ord , ord)  ; 

7.  Read  in  the  data  from  the  file  and  store  in  a  1XN  matrix 
EGM96=double (dlmread ( 1 EGM96 . dat ; ) ) ; 

7.  Used  for  checking  my  values  to  the  original 
7.  fid  =  fopen(’  EGM96_grob.txt  ’  ,  ’wt  ’  )  ; 
l 

7.  for  m  =  1  :  500 
7.  for  n  =  l  :  2 

7.  fprintf  (f  id  ,  '7.1  .  Of  ’  ,  EGM96  (m  ,  n)  )  ; 

7.  end 

7.  for  n  =  3  :  4 

7.  fprintf  (fid  ,  ’7.1  .  He  ’  ,  EGM96  (m  ,  n)  )  ; 

7.  end 

7.  fprintf(fid,’\n’); 

7.  end 
7. 

7.  fclose(fid); 

7«Define  mu  and  radius  of  the  earth 
mu  =  EGM96  (  1  , 1 )  ;  7,  m"3/sec~2 
R_Earth  =  EGM96  (2 , 1)  ;  7.  m 

7.  Do  a  for  loop  to  populate  the  C  and  S  matrices 
7.  with  the  data  from  the  file 
n  =  3  ; 

for  i = 1 : ord 

for  j  =1 :  i 

C(i,j)=EGM96(n,3)  ; 

S(i,j)=EGM96(n,4)  ; 
n=n+ 1 ; 

j=j  +i ; 

end 

i=i+l ; 

end 

7.  Unnormalize  the  C  and  S  coefficients 
for  k=0 : ord - 1 
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for  1 =0 : k 

if  1==0 

eps  =  1 ; 

else 

eps  =2 ; 

56  end 

fact=(eps*(2*k+l) ^factorial (k-1) ) / factorial (k+1) 
C(k+1 , l+l)=double (sqrt (fact)*C(k+l , 1+1) ) ; 

S(k  +  1 , l  +  l)=double (sqrt (fact)*S(k+l , 1  +  1) )  ; 

end 

61  end 


7. 

7. 

7. 

66  7. 
V. 
1 
% 
% 

71  7. 

7. 

7. 

% 


Used  for  checking  my  values  to  the  original 
fid  =  fopen(’EGM96_grob.txt  ’  ,  ’ wt  ’)  ; 


for  n=0 : ord - 1 
f or  m  =  0 : n 

fprintf  (fid  ,  ’ "/ 1 . 0  f 
fprintf  (fid  ,  ’  7, 1 . 0  f 
fprintf  (fid,  ’70l.lle 
fprintf  (fid,  ’  7,1. lie 

end 

end 

f close (fid)  ; 


’  ,n)  ; 

’  ,m)  ; 

’ , C (n+1 , m+1 ) ) ; 

\n ’ , S (n+1 , m+1) ) ; 


76  end 
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B.2  Inertial  Angular  Momentum,  Matlab  Code 


Listing  B.3:  Inertial  Angular  Momentum 

(appendix2  /  angularmom .  m) 

function  angularmom 

7.  This  function  reads  in  time,  position  and  velocity  data 
7,  from  a  text  file,  calculates  the  Inertial  Angular  Momentum 
4  7,  and  outputs  the  z-component  to  a  text  file. 

7«  Put  the  data  from  eci.txt  into  new  matrices  time,  x  and  xdot 
clc 

In  =  dlmr ead (’I:\My  Document  s\Thesis\EarthSat\Data\Circular30LongTBP. 
\ eci . txt  ’  )  ; 

9  limit  =  length  (  In)  ; 

a  =  l ; 
b  =  l ; 

while  b  <  limit 
14  time (a)  =  In (b)  ; 

t ime=double ( t ime ) ; 
a=a+ 1 ; 
b=b+3 ; 

end 

19 

c  =  l ; 

d  =  2; 

while  d  <  limit 

x (c , 1) =In (d  ,  1)  ; 

24  x (c ,2) =In(d  ,2)  ; 

x  (c  ,  3)  =In (d  ,  3)  ; 
x  =  double  (x)  ; 
c  =  c  +  l  ; 
d=d+3 ; 

29  end 

e  =  l ; 
f  =3 ; 

while  f  <  limit+1 
34  xdot(e,l)=In(f,l); 

xdot(e,2)=In(f,2)  ; 
xdot(e,3)=In(f,3)  ; 
xdot=double (xdot ) ; 
e  =  e  + 1 ; 

39  f=f+3; 

end 

limit2=limit/3; 

fid  =  f open( 1 1 :\My  Document s \ Thes i s \ EarthSat \Data\ ..  . 
Circular30LongTBP\Lz_grob . txt ’ ,  ’wt ’ ) ; 

44  fprintf (fid  ,  ’ Time  (TU)  Z  component  of  angular  momentum\n ’ )  ; 

fprintf  (fid,  ’ - \n  ’  )  ; 

for  j  =1 : 1 imit2 


xl=[x(j,l)  x ( j  ,2)  x ( j  ,  3) ]  ; 
xdot 1 = [xdot ( j , 1 )  xdot(j,2)  xdot(j,3)]; 
49  H=cross(xl,xdotl); 

fprintf(fid,  1  7,1 .6  f  ’  ,time(j))  ; 

fprintf  (f  id  ,  ’  7.0 . 15  f  \n  ’  ,  H  (3)  ) 

end 

fclose (fid)  ; 

54  f printf ( ’ D0NE\n  ’  )  ; 
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B.3  Change  in  Ascending  Node  Matlab  Code 

Listing  B.4:  Nodal  Regression 

(appendix2  /  nodalregress.m) 

function  nodalregress 

7„  Put  the  data  into  new  matrices  time,  x  and  xdot 
clc 

In  =  dlmr ead (’I:\My  Document s\Thesis\EarthSat\Data\Circular30LongTBP .  .  . 

\ eci . txt  1  )  ; 
limit=length ( In) ; 

a  =  l ; 
b  =  1 ; 

while  b  <  limit 

time (a) =In (b) ; 
time=double (time) ; 
a=a+ 1 ; 
b=b+3 ; 

end 

c  =  l  ; 
d  =  2; 

while  d  <  limit 

x (c , 1) =In (d  ,  1)  ; 
x (c , 2) =In (d  ,  2)  ; 
x(c,3)=In(d,3)  ; 
x  =  double (x)  ; 
c=c+l ; 
d=d+3 ; 

end 

e  =  l  ; 
f  =3 ; 

while  f  <  limit+1 

xdot(e,l)=In(f,l)  ; 
xdot(e,2)=In(f,2)  ; 
xdot(e,3)=In(f,3) ; 
xdot=double (xdot ) ; 
e  =  e  + 1 ; 
f=f +3; 

end 

limit2=limit/3; 

fid  =  fopen(’I:\My  Document s \ The s i s \ EarthSat \Data\ ..  . 

Circular30LongTBP\Nodal_Regress_grob . txt ’ ,  1  wt ’ ) ; 

for  j  =1 : limit2 

xl  =  [x(j,l)  x ( j  ,2)  x ( j  ,  3) ]  ; 

xdot 1 = [xdot ( j , 1 )  xdot(j,2)  xdot(j,3)]; 

H=cross (xl , xdot 1 ) ; 

N=[-H(2)  H  (  1 )  0]; 
if  N(2)  <  0 

omega=2*pi  -  acos (N (1) /norm (N) ) ; 
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else 

omega=acos (N (1) / norm (N) ) ; 
end 

if  omega  <  0.0 

omega  =  2*pi  -  omega; 

end 

fprintf(fid,  ’  °/« 1 . 6  f  ’  ,time(j))  ; 
fprintf(fid,’  "/oO.lBfXn1, omega) 

end 

fclose (fid)  ; 
fprintf ( 1  DONE \n  ’  )  ; 
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B.4  Least  Squares  Matlab  Code 


Listing  B.5:  Least  Squares 

(appendix2 /leastsquares.m) 

function  leastsquares 

2  l  This  program  estimates  data  using  a  simplified 
%  dynamics  model,  creates  an  initial  estimate 
7„  using  linear  least  squares  and  then  updates 
°/0  the  estimate  using  non-linear  least  squares 

7  clc 

t  =0 : 60 : 600 ; 

7,  asssume  accuracy  to  meters 

sigma=.l;  %  Accuracy  of  the  sensor  and  noise  level 

12 

7.  Create  the  simulated  data 

•/.  Define  the  constants  and  ICs 
R0=  [1 . 05*6378137  ;  0  ;  0]  ; 

V0 =[0;6538. 656905;3600. 9857675]; 

17  g=-9.81; 

7,Run  a  for  loop  to  calculate  position 
for  i=l : length (t) 

R ( :  ,i)=R0  +  V0*t(:  , i )  +  [0 ; 0 ;  ,5*g*t(:  , i ) ~ 2]  ; 

22  end 

%  Add  random  noise  to  the  data 
P=rand (3 , length (t ) ) ; 
z  =  R  +  P*  sigma ; 

27 

y.  Run  linear  least  squares  to  find  an  initial  estimate 
“/,  Define  constants  and  initiate  sums 
Q=[sigma~2*eye (3 ,3) ]  ; 

H=[eye(3,3)  , zeros  (3,3)]  ; 

32  TQIT_sum=0 ; 

TQIz_sum=0 ; 

•/.  Run  a  for  loop  to  complete  the  two  summations 
for  i=l : length (t) 

STM=[eye(3,3) ,t(: ,i)*eye(3,3) ; zeros (3,3) ,eye(3,3)] ; 
37  T=H*STM ; 

TQIT=T ’ *inv(Q) *T ; 

TQIz=T ’ * inv ( Q ) *z  (  :  , i)  ; 

TQIT_sum=TQIT_sum+TQIT ; 

TQIz_sum=TQIz_sum+TQIz ; 

42  end 

•/.  Calculate  the  covariance 
P_x=inv (TQIT_sum) ; 

•/.  Calcalute  the  initial  state 
state_init=P_x*TQIz_sum ; 

47 

•/.  Run  non-linear  least  squares  to  caculate  the  accuracy 
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7„  Initiate  the  sums 
TQIT_sum2=0 ; 

TQIr_sum=0 ; 

52  */.  Run  a  for  loop  to  calculate  the  sums 
for  j =1 : length (t) 

STM=[eye(3,3) ,t(: ,j)*eye(3,3) ; zeros (3,3) ,eye(3,3)] ; 

"/«  Find  the  state  and  add  a  particular  solution 
state=STM*state_init+[0;0;0.5*g*(t(: ,j))~2;0;0;g*t(: ,j)] ; 
57  1  Treat  as  a  linear  problem 

G=[eye(3,3) , zeros (3,3)] *st ate; 

H2=[eye(3,3)  ,  zeros (3,3)]  ; 
r=z ( : , j ) -G ; 

T2  =  H2  *  STM  ; 

62  TQIT=T2 ’*inv(Q)*T2; 

TQIr  =  T2 ’ * inv (Q) *r  ; 

TQIT_sum2=TQIT_sum2+TQIT ; 

TQIr_sum=TQIr_sum+TQIr ; 

end 


67  °/,C  Calculate  the  covariance 
P_dx= inv ( TQ IT_ sum2 ) ; 

%  Calculate  the  initial  state  variance 
dstate_init  =  P_dx*TQIr_sum ; 

"/,  Calculate  the  updated  initial  state 
72  state_init=state_init+dstate_init ; 

•/.  Calculate  the  initial  state  error 

R_diff  =  [abs (RO (1) -state_init  (1)  )  ; abs (RO (2) -state_init  (2) )  ; abs (RO . . 
(3) -state_init (3) )] ; 

V_diff  =  [abs (VO (1) -state_init  (4)  )  ;  abs (VO (2) -state_init (5) )  ; abs (VO . . 

(3) -state_init (6) )] ; 

7,  Output  the  data  to  the  screen 

77  fprintf (’ Position  error  is  : \nx-direction  °/« 1  .  1 2  e \ny  -  dir  e  ct  i  on  °/01.12 
e\nz -direction  °/„l  .  12e\n ’  ,R_diff  (1)  ,R_diff  (2)  ,R_diff  (3)  )  ; 
fprintf  (  1  Velocity  error  is  : \nx-direction  °/0l  .  12e\ny-direction  °/01.12 
e\nz-direction  %1 . 12e\n’ ,V_diff (1) ,V_diff (2) ,V_diff (3) ) ; 
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